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Detection  of  changes  in  longitudinal  data  is  very  important  in  medicine  and 
industry.  However,  this  process  can  be  influenced  by  outliers,  leading  to  incorrect 
conclusions  regarding  the  changes.  The  purposes  of  this  research  were  to  identify 
a scheme  to  reduce  the  impact  of  outliers,  to  estimate  the  parameters  necessary  to 
implement  the  scheme,  and  to  investigate  its  effect  on  change  detection. 

To  reduce  the  impact  of  outliers,  an  option-3  scheme  is  proposed.  An  ad  hoc 
version  of  this  scheme  has  been  used  by  some  clinical  researchers  without  theoretical 
justification.  In  order  to  determine  the  properties  of  outliers,  the  errors  are  assumed 
to  come  from  a mixture  of  two  normal  distributions,  one  with  expected  small  vari- 
ation and  the  other  with  unusually  large  variation.  The  option-3  scheme  can  be 
justified  as  a robust  scheme  against  nonnormality  of  a large  variation  error  distri- 
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bution.  There  is  considerable  sample  size  gain  in  using  this  scheme  over  the  usual 
fixed  sample  schemes.  Since  the  small  error  variance  in  the  mixture  distribution 
was  crucial  to  the  option-3  scheme,  a new  estimation  method  for  the  variances  and 
proportions  in  a mixture  of  two  normal  distributions  is  developed.  When  sample 
size  is  small,  the  new  estimates  are  more  accurate  than  those  based  on  the  maximum 
likelihood  procedure. 

Many  methods  for  detecting  changes  in  longitudinal  data  were  studied  under 
the  most  commonly  used  periodontal  disease  progression  models,  gradual  and  burst 
models.  Under  these  two  models,  the  regression  method,  the  single  change  point 
method,  and  a mixture  of  the  two  methods  were  found  to  be  more  sensitive  in 
identifying  changes  than  the  other  methods.  The  option-3  scheme  has  been  applied 
to  these  three  change-detection  methods,  and  it  has  been  found  that  the  option-3 
scheme  detects  changes  more  quickly  than  the  fixed  two-sample  scheme. 
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CHAPTER  1 
INTRODUCTION 


Detection  of  changes  in  longitudinal  data  is  important  for  testing  hypotheses  on 
growth  patterns.  One  of  the  fields  involved  in  detecting  changes  in  longitudinal  data 
is  found  in  periodontal  disease  research.  A current  “gold  standard”  to  determine 
periodontal  disease  activity  is  the  change  in  attachment  level  observed  by  comparing 
the  probe  readings  during  a patient’s  regular  visits  to  a clinic.  However,  these  probe 
readings  may  not  always  be  accurate.  There  is  literature  in  recent  periodontal  jour- 
nals on  the  accuracy  of  the  attachment  measurements  (e.g.,  Watts,  1987;  Janssen  et 
al.,  1988;  and  Karim  et  al.,  1990).  It  has  been  shown  that  even  under  very  carefully 
controlled  clinic  conditions,  outliers  cannot  be  eliminated  (see  Goodson,  1986;  Gun- 
solley  and  Best,  1988;  and  Yang  et  al.,  1992).  Figure  1.1  gives  a typical  histogram 
of  measurement  errors.  It  is  constructed  based  on  samples  of  differences  between 
two  measurements  of  the  same  site.  Thus  the  histogram  is  centered  at  zero.  In  most 
periodontal  studies,  six  sites  per  tooth  are  probed,  and  any  clinically  significant  loss 
in  attachment  at  any  one  of  the  sites  is  considered  unhealthy,  requiring  surgical  or 
other  types  of  treatment  (Grbic  et  al.,  1991).  Since  there  may  be  192  sites  in  a 
mouth,  depending  on  the  number  of  teeth,  a small  percent  of  outliers  can  easily 
produce  a false  alarm. 
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Outliers  have  been  discussed  since  the  17th  century.  Beckman  and  Cook  (1983) 
summarized  the  history,  and  two  books  by  Barnett  and  Lewis  (1984)  and  Hawkins 
(1980)  contained  many  examples  and  techniques  about  outliers.  There  are  two 
mainstream  practices  on  manipulating  outliers,  one  of  which  is  to  remove  outliers 
which  require  filtering.  The  other  is  to  keep  them  as  a genuine  reflection  of  the  data. 
This  dissertation  is  concentrated  on  the  former  case.  Anscombe  (1960),  Guttman 
and  Smith  (1969),  Guttman  (1973),  Marks  and  Rao  (1978)  and  others  approached 
rejecting  outliers  by  using  various  estimation  techniques.  Initiated  by  attachment 
level  (AL)  measurements  in  periodontal  clinics,  the  following  protocol  is  suggested 
as  a method  to  remove  outliers  and  to  get  a better  estimate  of  true  AL  change:  Two 
probings  are  taken  initially,  and  if  the  two  measurements  are  not  close  enough,  a 
third  measurement  will  be  taken.  Some  clinical  researchers  have  already  used  this 
proposed  scheme  or  a similar  one.  However,  there  is  a lack  of  consensus  on  how 
close  is  close  enough.  Osborn  et  al.  (1990)  chose  the  threshold  as  1mm,  Best  et 
al.  (1990)  chose  3mm,  while  Marks  et  al.  (1991)  chose  1.1mm  for  their  clinical 
periodontal  studies. 

The  main  purpose  of  this  dissertation  is  to  study  the  optimal  choice  of  this 
threshold,  to  estimate  the  parameters  necessary  to  implement  this  threshold,  and 
to  investigate  its  effect  on  change  detection.  A brief  description  of  the  main  idea 
presented  in  this  dissertation  is  as  follows.  Let  x\  and  x2  be  the  first  two  measure- 
ments of  a certain  quantity.  If  \x\  — x2|  < 8 for  some  given  6,  then  mean  estimate 
A = (xi  + m2)/2;  otherwise,  a third  measurement  will  be  taken  and  the  mean  is 
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estimated  by  the  average  of  the  closest  two.  This  new  rule  is  named  the  “option- 
3 scheme.”  To  determine  6,  mean  squared  error,  E(fi  — p)2,  is  used,  which  is  a 
function  of  error  distribution. 

To  describe  the  distributions  that  produce  outliers,  heavy-tailed  distributions 
were  considered  first  (Daniell,  1920;  Fisher,  1922;  Student,  1927).  Later,  in  the 
20th  century,  mixture  models  were  suggested  by  Ogrodnikoff  (1928),  Jeffreys  (1932), 
Dixon  (1950),  and  Grubbs  (1950).  The  most  common  mixture  assumption  is  that 
the  error  follows  a mixture  of  two  normal  distributions:  one  with  variance  a2  and 
the  other  with  variance  k2cr2,  where  k2  > 1.  Due  to  the  scarcity  of  outliers,  its  dis- 
tribution may  not  be  easy  to  identify.  Hence,  the  robustness  against  nonnormality 
of  the  outlier  distribution  is  also  investigated. 

The  values  for  optimal  6 depend  on  the  error  distribution,  especially  the  variance 
a2  and  its  proportion  p in  the  mixture.  Thus,  the  parameters  cr2  and  p should  be 
estimated  as  accurately  as  possible.  The  maximum  likelihood  estimation  method  to 
determine  variances  and  the  proportion  in  the  mixture  of  two  normal  distributions 
has  been  discussed  previously  by  Behboodian  (1970),  Hosmer  (1973a),  Day  (1969), 
and  others.  However,  these  estimates  do  not  give  good  estimates  when  the  sample 
size  is  small.  The  reason  is  that  the  two  components  of  the  mixture  distributions 
are  not  well  separated  (Hosmer,  1973b).  A new  estimation  method  is  developed  for 
more  accurate  estimation  of  a2  and  p when  the  sample  size  is  small. 

Several  change  detection  methods  for  longitudinal  data  are  being  used  in  pe- 
riodontal disease  studies.  Among  them  are  running  median,  regression,  tolerance 
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methods  (Haffajee  et  al.,  1983),  and  cusum  (Aeppli  and  Pihlstrom,  1989).  Aep- 
pli  and  Pihlstrom  (1989)  claimed  cusum  and  running  median  methods  are  better 
than  regression,  while  Yang  et  al.  (1991)  claimed  that  regression  is  better  in  most 
cases.  The  discrepancy  between  the  two  conclusions  is  due  to  differences  in  the 
assumption  of  error  distribution,  i.e.,  Aeppli  and  Pihlstrom  (1989)  did  not  use  the 
same  assumptions  for  all  three  methods;  measurement  error  distribution  is  known 
for  running  median  and  cusum,  but  not  for  regression.  On  the  contrary,  Yang  et 
al.  (1991)  used  the  same  assumptions,  with  measurement  error  distribution  being 
known  for  all  methods.  Another  difference  between  the  two  articles  is  the  error  dis- 
tribution; logistic  error  distribution  was  used  in  Aeppli  and  Pihlstrom  (1989),  while 
normal  error  distribution  was  used  in  Yang  et  al.  (1991).  These  discrepancies  are 
reconsidered  in  this  study.  It  is  found  that  the  regression  method  identifies  changes 
more  quickly  than  the  two  other  methods,  both  for  logistic  and  normal  error  dis- 
tributions. Four  new  methods  for  detecting  changes  are  introduced  and  compared 
with  the  regression  method  under  the  normal  error  structure.  The  change  detection 
methods  with  option-3  scheme  is  also  studied  in  change  identification. 

Chapter  2 presents  the  statistical  model  for  outliers’  distribution.  Mathematical 
formulation  of  the  option-3  scheme  is  shown,  and  some  advantages  of  the  option-3 
scheme  are  discussed.  The  result  obtained  by  applying  option-3  to  real  data  is  pre- 
sented. Chapter  3 discusses  a new  estimation  method  of  parameters  for  the  option-3 
scheme.  A comparison  of  the  new  estimation  and  standard  maximum  likelihood  es- 
timation is  given.  Chapter  4 describes  currently  used  methods  for  change  detection 
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in  longitudinal  data.  They  are  compared  with  some  newly  developed  methods. 
Chapter  5 discusses  the  interaction  between  change  detection  and  option-3  mea- 
surement schemes.  Chapter  6 summerizes  the  results  and  provides  a brief  discussion 
of  future  research  issues. 
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Figure  1.1:  A histogram  of  probing  measurement  errors 
ment  error) 


(frequences  vs  measure- 


CHAPTER  2 

OUTLIER  REDUCTION  BY  AN  OPTION-3  MEASUREMENT  SCHEME 

2.1  Introduction 

As  we  reviewed  in  Chapter  1 , outliers  can  limit  the  ability  of  our  conclusion  and 
decision.  This  chapter  deals  with  identifying  and  handling  outliers.  To  accomplish 
this  goal,  we  consider  the  following  scheme.  If  the  first  two  measurements  are 
close  enough,  analysis  is  carried  out  with  those  two  measurements,  but  if  the  first 
two  measurements  are  not  close  enough,  a third  measurement  will  be  taken.  More 
specifically,  if  the  difference  of  the  first  two  measurements  is  smaller  than  a given 
value,  say  8,  then  the  average  of  the  first  two  measurements  is  taken  as  an  estimate 
of  the  true  value.  But  if  the  difference  is  larger  than  8,  then  a third  measurement 
is  taken  and  the  average  of  the  closet  two  values  among  the  three  measurements  is 
regarded  as  the  estimate.  Let  us  define  this  rule  as  the  “option-3  scheme.”  The 
main  issue  in  the  option-3  scheme  is  how  to  determine  the  threshold,  8.  Various 
values  of  8 have  been  chosen  in  the  past  due  to  lack  of  theoretical  support.  Osborn 
et  al.  (1990)  chose  the  threshold  as  1mm,  Best  et  al.  (1990)  chose  3 mm,  while 
Marks  et  al.  (1991)  chose  1.1  mm  in  periodontal  disease  studies.  Hence  the  main 
point  of  this  chapter  is  to  find  the  optimal  8. 

For  studying  the  property  of  outliers,  error  distribution  is  assumed  to  be  a mix- 
ture of  two  normal  distributions.  Mixture  models  have  been  suggested  for  outliers 
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by  Daniell  (1920),  Ogrodnikoff  (1928),  Jeffreys  (1932),  Dixon  (1950),  and  Grubbs 
(1950).  To  find  the  optimal  8 , the  mean  squared  error  is  used  to  measure  the  dis- 
crepancy between  the  true  value  and  the  estimated  value.  The  numerical  solution 
of  8 is  obtained  for  the  mixture  of  two  normal  distributions.  Due  to  the  difficulty 
in  analyzing  distributions  other  than  normal,  we  use  a simulation  study  to  examine 
the  robustness  of  the  option-3  scheme  against  the  nonnormality  of  outliers. 


2.2  Model  and  Method 

Let  N(p,  a2)  denote  the  normal  distribution  with  mean  p and  variance  a2,1  and 
Nk,q(p,cr2)  represent  the  mixture  of  two  normal  populations  with  a proportion  p 
from  N(p,cr2)  and  q(=  1 — p)  from  N(p,k2a2).  Let  x\,x2,X3  be  independent  and 
identically  distributed  (i.i.d.)  random  variables  from  Nkiq(p,cr2). 

Suppose  <7  = 1.  Then  the  model  can  be  expressed  as 

Xi  = p + ei  i-  1,2,3 


where  e*  is  distributed  as  a mixture  of  iV(0, 1)  and  N( 0,  k2)  with  probability  density 
function 


/(e)  = 


P 


-.e  2 + 


— g_ 

e 2^ 


y/2ir  \/2irk 

The  option-3  estimate  of  p,  denoted  by  fi,  is  given  by 


(2.2.1) 


*1+S2 

2 


sum  of  closest  two 
2 


if  I®!  — x2\  <8 


if  |xx  — x2\  > 8 . 


(2.2.2) 


^(0, 1)  distribution  is  called  standard  normal  distribution. 
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Equation  (2.2.2)  depends  on  the  threshold  value  6 , which  will  be  chosen  so  that  it 
minimizes  E(£l  — /x)2  (=  var(/x),  since  E(/x)  = /x  can  be  proven  by  symmetry).  The 
mean  squared  error  of  /x  in  (2.2.2)  is 


-£(£  ~ vY  - E [(£  - /x)2/|Xl_X2|<5  + E [(/x  - m)2/|Xi_*2|>5 


(2.2.3) 


where  I a is  the  indicator  function  which  is  1 if  A is  satisfied  and  0 otherwise.  Now 
for  the  general  a case,  we  let  2/i , 3/2, 3/3  be  random  samples  from  iVjt>q(/x,  a2).  Let 

if  |3/i  ~ 2/2I  < S' 
if  \yi  -3/2I  > S'  . 


VI+V2 

2 


(2.2.4) 


sum  of  closest  two 
2 


Then  it  can  be  shown  easily  that 


E(/X*  - /x)2  — C2E  [(/X  - /x)2/|Xl_Xj \<S'/tr  + [(£  - /i)2-f|xi-x2|>5'/<7  • (2.2.5) 

From  equation  (2.2.5),  it  is  found  that  the  threshold  is  aS  when  a is  not  equal  to 
one.  Thus,  without  loss  of  generality,  let  a = 1.  Analytic  computation  of  E(£l  — fj,)2 
and  simulation  are  explained  as  follows. 

The  first  part  in  (2.2.3)  can  be  represented  as 

E [(£  - /i)2f|x1-x2| <s]  = JJ  ^—--~2  - f(x1)f(x2)dx1  dx 2. 

Let  Zi  = Xi  — fi,  i = 1, 2,  then 

E [(£  A1)  -f|xi-x2|<5 


= ttjU  K.  (H^)2  (2-2.6) 
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where 


Let 


a,  = i 


i = 1 


— i — 2 
2 1 ~ z 


p -(-i-M)2  9 .-ta-e£ 

/(li)  = vr  1 +Vl5e  ’ 


5(z0  = 


„ _2  _ _ .2 

p zil  , q _i 

* 1 — — - — 2 k* 


, — -e  2 -\ — -'  ev?  , for  i = 1,2. 

V2^  y/2*k 


<H1(a,b,a,/3)  = j*  e-^-^dx 

= ^ [«(V^(6  - m - $(Vm«  - /?)] , 


$2(a,  6,  a,  0)  = [b  xe-a(-x-W dx 
J a 


2 a 


,-a(a-P)J-a(6-0)2 


+ /3'Jf1(a,&,a,/3), 


= [b  x2e~a(x-0)7dx 

Ja 

= ^ [(°  + P)e-a(a-0)2  ~(b  + P)e~a{>-^]  + 

^^x(a  - P,  b - p,  a,  0)  + P2Vi(a,  b,  a,/3) , 

where  $ is  the  cumulative  standard  normal  distribution. 

Let  t = Z\  + Z2  and  u = z\  — z2,  then  the  ranges  of  new  variables,  t and  it,  are 
|u|  < 8 and  — oo  < t < oo. 


2 , 2 
aiZ1  + ajZ2 


And 


11 


ft-\-u\2  ft  — u\'4 

= «(— ) +«<(—) 


1 


— t s (a*  + aj)  ( t + 

4 I \ ai  + aj 

— —hiCt-\-  u\  + h3u2 


(di  — dj)u\2  ( di  — dj)2U2 


(di  4-  ay) 


+ (di  + dj)u2 


where  hi  = di  + ay,  h2  = a*  — ay,  h3  = (ayay)//ii.  So  a typical  term  in  (2.2.6)  can 


be  represented  as 


ff  le-^4)dzidz2 

J J\zi—  zi\<6  V 2 / 7T 


= £El[‘  r 

07T  J-SJ-oo 


e~hsU  du 


/didj 


Ay/jrh^2 

From  (2.2.6)  and  (2.2.7),  the  first  part  of  (2.2.3)  equals 
[(A  — A4)  -f|*i-x2|<i] 


|2/i1^rx(— 6, 6,  h3,  0)  -f  h\  ^3(— 8,  S,  h3, 0)j  . (2.2.7) 


2 2 


= £ £ 77=7iW  6,  h3, 0)  + h\  V3(-S,  6 , fca,  0)}  . (2.2.8) 

»=i  y=i  4A/7rh1/ 

The  second  part  of  (2.2.3)  can  also  be  expanded  as  follows. 

= / fl^>s  (r--p- - „)2  /(xO/txO/^Mx,  <fa2  <fe3 
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= /// 


|xi-x2|>5  i 2 

|xi -xj |<mtn(|xi -x* |,  |x2-x3|) 


+ 

+ 


ill 

III 


|xi  -x2 |>£  V 2 

|xi-xs|<m»n(|xi-x2|,  |x2-x3|) 


lxl-*2l>5 

x2-x3|<mtn(|xi-x2|,  |xi.-x3|) 


^Zi_+_X2  _ f(x1)f(x2)f(x3)dx1  dx2  dx3 

^Si_+X3  _ ^ f(Xl)f(x2)f(x3)dxi  dx2  dx3 
('X2_+_X3  _ ^ f(x1)f(x2)f(x3)dx1  dx2  dx3 


- Ill 


\zl-z2\>S  \ 2 

|*l-x2|<m»'n(|zi-x3|,  |z2— z3  |) 


+ 

+ 


III 

III 


(~ -y^2)  g(zi)g(z2)g(z3)dzi  dz2  dz3 

(Zl  ^ Z3)  5(21)5(22)5(23)^21  dz2  dz3 


|zi-z2|>i  \ 2 

|zi-zj|<m«n(|zl-z2|,  |z2-z3|) 

|,1-*2|>«  (~2  2 ^ ) 5(2l)5(22)5(23)<i2i^2dz3 

|z2-z3|<mtn(|zl-z2|1  |zi-z3|) 


— Ai  + A2  + A3 

where  Ai,  A2,  and  A3  denote  the  first,  the  second,  and  the  third  term,  respectively, 
and  f(xi)  and  5(2^)  are  defined  the  same  way  in  the  first  part  of  the  integration. 
Ai  can  be  written  as 


* “ E ///  h - W i^r) 2 

*tj<k  |zl-z2|<m»n(|zi-z3|,  |z2-z3|) 

where  a»,  a,j,  and  a*,  are  defined  the  same  way  in  the  first  part  of  the  integration. 
Let  v = Z\  + z2,  x — z\  — z2,  y = 21  + z2  — 2 z3.  Then  the  range  of  the  new  variables, 
v,  x,  and  y,  is  v U R where  v = (—00,  +00),  and 


(Rl) 

y > 36, 

< x < -8 

(R2) 

y > 38, 

8 < x < | 

(R3) 

y < -36, 

8<x<~l 

(R4) 

y < -38, 

\<x<-8. 
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Also 


axz\  + a,jzl  + akz\ 


( v + ®\ 2 (v  — x\2  (v  — y\'‘ 

*(— ) +■*(—)  +a‘(_21) 


= 7 (<■<  + + .,)  (v  + ^ ~ °j)l 

4 1 \ dx  4-  a,j  -\-  dk  ) a;  4-  aj  4-  dk 


(ai  + a_,)x2  + afcy2} 


= l|9(v+(ai  +7(*  + Avf  + OTJ 


where  /?  = a;  + a^-  + a*;, 
4aiaj  + (a<  + a.,)afc 


7 = 


4 P 


(a»  aj)ak 

> A = 7“ — ~ r~~T~»  S' 


diCLjd  k 


Adidj  + (di  + dj)dk  ’ 4didj  4-  (ai  + o.j)dk 


So  a typical  term  in  A\  is 


/// 


^a“l+‘“i+a-i)^dzIdz3 

|zx-z2|<m»n(|zi— zs|,  |z2-zs|) 


= jx  fff  v2e  Kv+(  ‘ **"  **)  e~^x+Ay)i e~9y2 dvdxdy 
16  J J JvUR 


= J JR  (2/3  + (a»  - aj)2x2  - 2{di  - dj)dkyx  4-  a£y2)  e 7(x+Ay)2  e ffy2  dxdy 


A 

8/35/2 


too  ty/3 

H dxdy  4-  / / U dxdy 

-v/3  >/« 

/•-34  f-y/3  t-36  t-S  1 

+ / / tl  dxdy  + / i dxdy  \ 

J — OO  J6  J — OO  Jy/3 


(2.2.9) 


The  first  term  inside  the  curled  brackets  above  can  be  represented  as 


f j D dxdy  = f /(— y/3,  —6, 7,  — Ay)e  sy2dy 

./3£  J-y/3  J3S 
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where 


/(-y/3,-5,7,-Ay)  = (0^  - aj)2^3(-y/3,-^7,-Ay) 

- 2 (ai  - aJ)afcj/^2(-y/3,  -5,7,  - Ay) 
+ (a£y2  + 2/?)'Pi(— y/3,  -5,7,  -Ay)  . 


Other  terms  can  be  represented  in  a similar  manner. 
Thus, 


A\ 


„ /aiaTaZ  ( , 

12  87T/3V2  {Jx  (/(-y/3,-5,7,-Ay)  + /(5,y/3,7,-Ay))e"av  dy 
+ J ^ (/(^-y/3,7,-Ay)  + /(y/3,-5,7,-Ay))e“3y2dy|.  (2.2.10) 


Now  consider  A2,  we  have, 


- ±w 

*,J,k  I. 


1*1— *2  I >6 

|*l-2j|<mm(|zi-*2|,  |*2  *3 1) 


(fL^)  2 )dzxdzidz* 


Let  t — z\  — z3,  v = z\  + z3,  u = z3  — z2.  Then  the  range  of  new  variables,  t , v, 
and  it,  is  v U S where  v = (— oo,  -foo)  and  S is 


(SI) 

| < u < 25, 

—u  + 5 < i < u 

(S2) 

u > 25, 

-!<*<“ 

(S3) 

-25  < u < - J, 

u < t < —u  — 5 

(S4) 

u < —25, 

«<*<-?• 

Also, 


a.iz\  + ajzf  + akzl 


ft  + v\2  ( v — t — 2u  \ 2 fv  — t\ 

= “>{—)  +“i(  j J +at(—J 
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1 , . , X / , (<Xi  — — afc)t  - 2 (LjU 

T(<*i  + Oj  + afc)  \ v + 3— f — 

4 ( aj  + ay  -}-  a*. 

aj(ay  + afc) 


t + 


\ 2 2 
a,u  \ ayafcU 


+ 


(aj  + ay  + afc)  \ ay  + a*  J aj  + a* 

(3  ( nt  — 2 a,u\2  . , s,  , , 

= ~ ( v H ^ J + a(4  + ^u)2  + ku2 


where  /3  = aj  + ay  + at,  n = ax  — a}  — a*, 


a = 

ai(aj  + ajt) 

h Gj 

J ayafc 

/?  ’ 

aj  + O'k 

ay  + ajfe 

So  the  typical  term  in 

A2  is 

///  . 

l*i  — *2 1>* 

(Z\  + Z3 > 

^ e-(^l+^zl+^l)dzldz2dz: 

1*1  — *3  |<rr»*r»(|*i  — *2 1,  \z2-za  |) 


= \ [[[  v2e  *(v+^~)  e-«(«+M 3e-ku*dvdtdu 
8 J J Jvus 

= / j>  j (nt  ~ 2ajU)2) 

H 


4/JS/2 


6 f2J — v-f*£ 


HD  dtdu  + J J DD  dtdu 

r—S/2  f—u—S  [—26  r—u/2  ~) 

^ J 26  J ^ J J ^ | ' 


(2.2.11) 


The  term  inside  the  curled  brackets  in  (2.2.11)  can  be  represented  as 


r2 6 [u  r26 

/ / D $ dtdu  = / p(—  u -f  6,  u.a,  — hu) 

J6/2J-U+6  J6!  2 ' 


Ar26 

6/2  ‘ 


e~ku'dy 


where 


p(—u  + 8,u,a,—hu)  = n2^3(— u + 8,  u,  a,  — hu) 


— 4 ajnu^2(—u  + 6,  u , a,  — /m) 

+ (4 a2jU2  + 2(3)'£i(—u  + 6 , u,  a,  — hu)  . 
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Similarly  other  parts  of  the  above  integration  can  be  represented. 
Thus, 

+ f p(—u/2,u,a,—hu)ekxi7du 
J 26 

r26  k 2 

+ / p(u,  — u — £,  a,  —hu)e*u  du 
J-S/2 

+ J p(u,—u/2,a,—hu)eku2dv)j  . 


(2.2.12) 


Next,  consider  A3.  we  have, 

M = ±JH  , ,l2t 

|*2-za|<m»n(|z1-z2|1  |«i  — *3 1) 

Let  t = z2  — z3,  v = z2  + z3,  u = z3  — Z\.  Then  the  range  of  the  new  variables,  t,  v , 
and  u,  is  v U S where  v and  S are  defined  in  p.14. 

Also, 


2 . 2 , 2 
a{zx  + a,jZ2  + ckz3 


= a,  ^ 


v — t — 2 u\2 

—5—)  * 


which  is  the  same  corresponding  form  in  A2  except  that  a;  and  aj  are  interchanged. 
Thus  A3  can  be  computed  easily  by  interchanging  a*  and  aj  in  A2.  With  (2.2.8), 
(2.2.3)  is  now  computed  by  single  integrals. 

Due  to  the  complexity  of  the  computation,  a solution  has  to  be  found  by  nu- 
merical integration.  For  selected  values  of  q (0.01,  0.025,  0.05,  0.1),  and  k (3,  4, 
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5),  Figure  2.1  shows  the  variances  of  [l  as  functions  of  8.  Figure  2.1  shows  that 
as  q increases,  the  optimal  8 decreases  for  a fixed  k value.  Also,  as  the  value  of 
k increases,  the  optimal  8 decreases  slightly  for  a fixed  value  of  q.  Variance  and 
optimal  6’s  are  shown  in  Table  2.1.  Also  the  optimal  8’s  are  plotted  in  Figure  2.2 
for  the  different  k and  q. 

Figure  2.1  and  Table  2.1  have  also  been  checked  by  simulation.  In  the  simula- 
tion study,  random  numbers  were  obtained  by  generating  U(0,1)2  r.v’s  using  Park 
and  Miller’s  (1988)  method  and  transformed  to  independent  N(0,1)  r.v’s  by  using 
the  Box-Muller  transformation  (Box  and  Muller,  1958).  A mixed  normal  popu- 
lation (2.2.1)  was  used  to  generate  data.  The  observations  were  assigned  to  the 
components  by  generating  another  U(0,1)  r.v.,  u.  If  u < p , then  the  observation 
was  assigned  to  a N(0,1),  otherwise  the  observation  was  assigned  to  a N(0,  k2). 
Variances  of  fi  with  the  corresponding  value  of  optimal  8 based  on  5000  simulations 
are  given  in  Table  2.2  for  various  combination  of  values  of  q and  k.  The  values 
agree  very  well  with  the  theoretical  ones  (Figure  2.1).  For  example,  when  k=3  and 
q=0.1(or  k=5  and  q=0.01),  the  simulated  optimal  value  of  8 is  around  3.0(or  4.0), 
very  close  to  the  theoretical  value. 


2U(0,1)  represents  uniform  distribution  of  interval  (0,  1) 


variance 
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K-3  K=4  K=5 


delta  delta  delta 


Figure  2.1:  Variance  of  mean  estimate  (2.2.2)  for  different  k and  q 
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Table  2.1: 


Minimum  variance  of  fi  and  corresponding  8 


k 

q 

8 

Var^p) 

3 

0.01 

4.4500 

0.524672 

0.025 

3.9500 

0.556232 

0.05 

3.5500 

0.606610 

0.075 

3.3000 

0.658318 

0.1 

3.1000 

0.713163 

4 

0.01 

4.3000 

0.526929 

0.025 

3.8000 

0.561301 

0.05 

3.4000 

0.620660 

0.075 

3.1500 

0.687690 

0.1 

3.0000 

0.764354 

5 

0.01 

4.2500 

0.526655 

0.025 

3.7500 

0.562736 

0.05 

3.4000 

0.631859 

0.075 

3.1500 

0.717237 

0.1 

3.0000 

0.820660 

Optimal  delta 
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Figure  2.2:  Optimal  value  of  8 for  different  k. 
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Table  2.2:  103x  Variance  of  jl  based  on  5000  simulations  of  samples  from  mixed 

normal  population 


k 

q 

O", 

II 

o 

o 

1.0 

2.0 

3.0 

4.0 

5.0 

6.0 

7.0 

8.0 

3 

0.01 

653 

641 

597 

562 

536 

538 

544 

548 

551 

0.05 

706 

695 

655 

630 

621 

646 

672 

688 

701 

0.1 

794 

785 

750 

734 

743 

791 

837 

868 

891 

4 

0.01 

652 

641 

596 

563 

538 

542 

552 

562 

570 

0.05 

714 

704 

665 

643 

643 

673 

719 

764 

802 

0.1 

838 

830 

796 

787 

812 

868 

952 

1033 

1105 

5 

0.01 

651 

640 

596 

562 

537 

543 

555 

567 

580 

0.05 

725 

715 

676 

655 

655 

694 

744 

803 

862 

0.1 

896 

888 

855 

846 

871 

943 

1034 

1141 

1248 
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2.3  Merits  of  Using  the  Option-3  Scheme 

This  section  describes  advantages  of  the  option-3  measurement  scheme.  The 
merits  discussed  here  are  an  advantage  in  sample  size,  robustness  to  outlier  distri- 
bution, and  closeness  to  normality. 


2.3.1  Sample  Size  Comparison 

An  obvious  question  arises  on  the  efficiency  of  the  estimate  (2.2.2).  The  average 
sample  size  of  this  estimate  is  greater  than  2.  Thus,  one  might  question  whether 
the  gain  in  accuracy  of  this  estimate  can  compensate  for  the  sample  size  increase. 
Let  X2, . . . , xn  be  i.i.d.  with  known  variance  cr2  and  let  the  variance  of  a mean 
estimate  fi  be  cr?.  If  the  mean  estimate  is  x,  then  <r?  = trj/n,  or  cr^/cr?  = n.  Let 
the  average  sample  size  used  to  construct  fi  be  n ^ and  define 


sample  size  benefit  of  fi  = — j — n £. 


(2.3.1) 


Obviously  a positive  value  of  equation  (2.3.1)  indicates  a gain,  and  a negative  value 
indicates  a loss  in  terms  of  sample  size.  Suppose  the  mean  estimate  is  x(=  J2x/n) 
for  fixed  sample  size  n,  then  a2~  = var(x ) = cr^/n.  Since  n#  is  n,  sample  size  benefit 
becomes  zero,  i.e.,  there  is  no  gain  or  loss  in  sample  size.  If  x is  based  on  two  (or 
three)  fixed  samples,  then  (2.3.1)  is  zero. 

Consider  the  estimate  of  mean  (2.2.2)  for  the  mixed  normal  population  (2.2.1), 
then 

a2  = var(x)  = p var  (x  ~ N(0, 1))  + q var  ( x ~ N( 0,  k2)j  = p -f  qk2. 
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For  instance,  for  k=3  and  q=0.05,  a2  — (0.95  + 0.05 *32)  = 1.4.  The  average  sample 
used  by  the  option-3  estimate  can  be  shown  to  be  equal  to 

rip  — 2Pr{|xi  — x2|  < <!>}  + 3Pr{|xx  — x2|  > £}.  (2.3.2) 


The  first  term  on  the  right  hand  side  of  equation  (2.3.2)  is 


Pr(|xi  - x2|  < 6) 


= Pr(— 8 < xi  — x2  < 6) 

/OO 

Pr(— 8 < xi  — x2  < 8\x2  = t ) Pr(x2  = t)dt 

-OO 

dxx  + —Lr-ef^)  dt 

h rk  ) V V2*k  ) 

J (probl[t  — 8,t  + 8])  (j)e~  + dt  , (2.3.3) 


-OO 

nt+S  i p q 

-S  Vv^tt  y/2irk 

l r°° 

2tT  J — oo 


where 


/■*+«  / _t2  q -t2  \ 

probl[t  — 8,t  + 8]  = J (^pe  2 + — e2*7 J 
The  second  part  of  the  equation  (2.3.2)  is 


dt 


Pr(|xi  — x2|  > 8) 

= f Pr(|xi  - x2|  > £|x2  = <)Pr(x2  = t)dt 

J — OO 

= f [Pr(xi  > t + 5)  + Pr(xx  < t - £)]  Pr(x2  = t)dt 

J — OO 


■ /:  \i 


t+S 


e 2 + 


■\Z2irk 


i ) dx\  + 


rt-S  / p 
J -oo  y y/2-K 


e 2 + 


y/2irk 


) dx\ 


( v =0.  q 

x ~7=e  2 + —/ ==r 
\ V27T  k 


eii?  1 dt 


/ -t2  q -t2\ 

= — / (probl[t  + 8,oo]  + probl[—oo,t  — 8])  [pe  2 + Te2^  ) 

2tt  J- oo  V k J 

Using  (2.3.3)  and  (2.3.4)  in  (2.3.2)  we  get 


dt 


(2.3.4) 


njx  = 2Pr{|xi  — x2|  < 6}  + 3Pr{|xi  — x2|  > 6} 
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{2  x probl[t  — 6,  t + 5]  + 3 x ( probl[t  + 6,  oo]  -f  probl[— oo,  t — 5])} 


X 


/ =£.  q =£\ 

lpe  1 + T ) 


dt 


(2.3.5) 


Figure  2.3  shows  the  sample  size  benefit  for  optimal  6 at  various  combination  of 
values  of  q and  k.  All  the  gains  are  positive  for  q,  i.e.,  the  option-3  scheme  is  more 
efficient  than  using  repeated  measurements  for  the  existence  of  outlier  measure- 
ments. 


2.3.2  Robustness  of  Outlier  Distribution 

Although,  as  one  can  see  from  Figure  1.1,  the  vast  majority  of  of  measurement 
errors  are  normal,  the  actual  distribution  of  the  outliers  is  difficult  to  identify.  To 
check  robustness  against  the  nonnormality  of  the  outlier  distribution,  two  distribu- 
tions are  studied:  a short  tail  uniform  distribution  and  a long  tail  double  exponential 
distribution.  The  error(e)  distribution  is  a mixture  of  a proportion  p from  N( 0, 1) 
and  the  remaining  proportion  q from  a double  exponential  or  uniform  distribution 
with  mean  0 and  variance  k2.  We  still  have  a2  — p + qk2. 

Figure  2.4  shows  the  sample  size  benefits  based  on  5000  simulated  samples  from 
the  mixtures  of  double  exponential  and  uniform  distributions  with  a N(0, 1)  distri- 
bution. With  the  exception  of  uniform  outlier  error  structure  with  small  q and  small 
k,  all  values  are  positive.  Outliers  with  a double  exponential  distribution  produce 
more  gain  than  normal  or  uniform  as  the  value  of  q increases. 


Sample  size  benefit 
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Figure  2.3:  Sample  size  benefit  for  optimal  8 based  on  option-3  rule 


sample  size  gain 
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Figure  2.4:  Sample  size  benefit  for  normal,  double  exponential,  and  uniform  error 
mixture  based  on  5000  simulations 
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2.3.3  Comparison  of  Kurtosis 


In  this  section,  we  try  to  see  whether  the  distribution  of  estimate  based  on  (2.2.2) 
is  closer  to  the  normality  assumption  than  that  of  usual  average  mean  estimate. 
Since  the  new  estimate  is  symmetric  about  /r,  the  skewness  is  the  same  as  for  a 
normal  distribution.  Kurtosis  is  used  to  test  the  deviation  from  normality. 

Let  Xi,X2, . . . ,x„  be  i.i.d.  with  mean  \i— 0 and  variance  a2.  Kurtosis  of  fi  is 
defined  E{fi)A/crA  — 3,  where  a2  = variance  of  fi.  For  the  usual  mean  estimate, 
xn  = (xi  + • • • + xn)/n,  the  kurtosis  value,  say  Kn , is  given  by 

K.-Z&- 3. 


Since 


E(K)  = 


Jr  (e  £(*?)  + Q E £(*?)£(*?) 

since  xjs  i.i.d.  and  E{x{ ) = 0 
= (n(Ki  + 3)cr4  + 3n(n  — 1)<74) 

Tb 

where  K\  is  kurtosis  for  one  observation, 
E(x  4) 

and  thus  K\  — r„,  — 3 


[E(x2)]2 


= —(Ki  + Zn) 


and 


E(x2)  = - , 
n 


kurtosis  of  xn  based  on  n observations  is 


E(K) 


Kn=  7^e^-3=  — 


[E(x2n)\ 


n 
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For  Nktq(0,cr2)  distribution, 


Kx 


E(*4) 

[E(x2)}2 


For  example,  k=5  and  p=0.9  gives  K\  = 
Define  that 


(p  * 3cr4  + q * 3fc4<74) 
(pa2  + qk2a2)2 

3Q  + qkA)  _ 

(. V + qk2)2 
13.45. 


reduction  gain  in  kurtosis  of  ft  = — rip, 


(2.3.6) 


where  n p is  the  average  sample  used  to  construct  ft  and  Kp  is  the  kurtosis  of  ft. 
Again,  a positive  value  of  the  equation  (2.3.6)  means  gain  and  a negative  value 
means  loss  in  terms  of  sample  size. 

In  Figure  2.5  the  reduction  gains  in  kurtosis  of  ft  based  on  5000  simulations  are 
plotted  for  various  value  of  k and  q.  Most  reduction  gain  in  kurtosis  are  positive, 
meaning  that  the  distribution  of  the  option-3  estimates  is  closer  to  normality  than 
that  of  the  conventional  fixed  sample  estimate  ( x ). 


2.4  An  Example 

In  periodontal  disease  investigations,  two  measurements  of  the  attachment  level 
(AL)  are  usually  taken  at  the  same  site  at  each  visit.  The  average  of  two  mea- 
surements is  used  as  an  estimate  of  the  true  AL  of  that  site  at  that  time.  Figure 
2.6  (a)  shows  a histogram  of  the  measurement  errors  of  probing.  The  measurement 
errors  were  estimated  by  the  difference  of  two  measurements  and  included  a small 
percentage  of  outliers.  Since  most  of  the  optimal  8 values  are  between  3.0  and 
4.0,  the  middle  value  3.5  is  selected  as  the  appropriate  8 value  for  this  application. 


Reduction  gain  in  kurtosis 


29 


q 


Figure  2.5:  Reduction  gain  in  kurtosis  for  optimal  8 based  on  5000  simulations 
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(a)  Original  data 


0.1mm 

(b)  Modified  data 


0.1mm 


Figure  2.6:  Comparison  of  original  data  and  modified  data 
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Thus,  the  true  threshold  is  3.5a  and  a is  estimated  across  all  sites  by  using  the 
difference  of  two  measurements  at  each  site,  which  will  be  discussed  in  Chapter  3. 
A histogram  of  the  modified  data  based  on  (2.2.2)  is  shown  in  Figure  2.6  (b),  which 
eliminated  all  the  outliers. 


2.5  Conclusion 

The  option-3  scheme  was  used  to  handle  outliers.  To  define  the  property  of 
the  outliers,  a mixture  of  two  normal  error  distributions  is  used.  Minimum  mean 
squared  error  criterion  is  used  to  measure  the  estimation  error.  Numerical  solution 
for  the  threshold  8 for  the  option-3  method  is  provided  in  Figure  2.1,  Table  2.1, 
and  Figure  2.2.  The  results  have  also  been  verified  by  simulation  (Table  2.2).  The 
results  show  that  the  optimal  8 value  is  usually  between  3.0  and  4.0  unless  the 
proportion  of  outliers  is  very  small,  generally  less  than  0.01.  Thus  8 = 3.5a  may  be 
safely  used  in  practice.  Section  2.3  shows  the  option-3  scheme  has  advantages  in 
sample  size,  in  robustness,  and  in  closeness  to  normality.  It  was  also  shown  based 
on  real  data  that  the  scheme  appropriately  removed  outliers.  The  option-3  scheme 
can  be  used  in  any  medical  science  and  industry  study  in  which  the  data  exhibit 
measurement  unstability,  but  require  accurate  estimate  of  a true  mean. 


CHAPTER  3 

A NEW  ESTIMATION  METHOD  OF  PARAMETERS 
3.1  Introduction 

In  Chapter  2,  the  determination  of  6 involved  two  variances  and  the  proportions 
in  a mixture  of  two  normal  populations.  The  smaller  variance  is  more  crucial  to  the 
threshold,  a6,  decision.  In  this  chapter,  estimation  of  variances  and  proportions  is 
discussed  for  the  mixture  of  two  normal  distributions. 

Traditional  estimation  of  parameters  in  the  mixture  of  two  normals  includes  five 
parameters;  the  two  means  /zi,  p.2,  the  two  variances  <7j,  a\,  and  the  proportion 
p.  However,  we  are  interested  only  in  the  measurement  error  variances.  Thus,  if 
we  take  the  difference  of  two  repeated  measures  that  have  the  same  mean,  then  we 
have  a mixture  of  normals  with  means  equal  to  zero. 

Maximum  likelihood  estimators  of  each  population’s  variance  and  proportion  for 
the  mixture  of  two  normal  populations  have  been  discussed  previously  by  Behboo- 
dian  (1970),  Hosmer  (1973a),  Day  (1969),  and  others.  However,  these  estimators  do 
not  provide  good  estimates  when  the  sample  size  is  small,  because  the  components 
of  the  two  distributions  are  not  well  separated  (Hosmer,  1973b).  A new  estimation 
for  the  variances  and  proportions  in  mixture  of  two  normal  distributions  is  proposed. 
An  example  based  on  real  data  is  used  to  demonstrate  the  estimation  of  the  pa- 
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rameters  of  error  distribution  and  the  reduction  in  variance  by  using  the  option-3 
scheme  with  the  estimated  parameters. 


3.2  Standard  Maximum  Likelihood  Estimation 


Let  the  mixture  of  two  normal  distributions  consist  of  N(0,  a j)  with  proportion 
p and  N(0,  of)  with  proportion  q(=  1 — p)  with  of  < of,  p » q. 

Let  Xi,  X2, . . . ,xn  be  independent  identically  distributed  random  variables  with 
density  function 

f(x)  = pfi(x)  + 9/a(x), 

where  0<p<l,<7  = l—  p, 

i =4 

/»(s)  = y— — - e 2lT»  , Z = l,2, 

V27 r<jj 

and  of  > 0,  i = 1,2. 

Then  log  likelihood  function  is 

L = lnft  /(xi)  = ]Cln  /(xj)  = Sln(p/ i(xi)  + C1  - P)Mxi)) 

3= 1 i=l  3=1 

Differentiation  with  respect  to  proportion  p and  ot2,  t = 1,2  gives 


dp 

dL 

dof 


h(xj)  h{xj) 


3=1 


/(X3) 


P 


/(*i)  \ 2\/2tto^ 


— V 1 — V 

i — g 


e + 


a;* 


%/2? Tdi 


§/(*<)!  2*lfl{X’)+2 

1 y-W  2 2\P/l(X3') 

’ 1 75T' 


(3.2.1) 


(3.2.2) 
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Similarly, 


Also, 


1 2 _ 2\(j;  P)f2(xj) 

do]  2oij£*  2)  fix,)  • 


ft  = 0 

dp 

^ = 0,  < = 1,2 


give  zero  value  of  the  left  hand  side  of  (3.2.1),  (3.2.2),  and  (3.2.3). 
Let  Pi  = p and  p2  = 9 and  define 


Wij=Pi 7TT  (i  = 1,2;  j = l,...,n). 

fixi) 


Multiplication  the  first  equation  in  (3.2.4)  by  p results  in 


Ph(xj)  ~ pMxj) 


3=1 


fiXi ) 


pfi(xj)  + (i  - p)h(xj)  - Mxj) 


3=1 


f(X3 ) 


Therefore 


Y A(*<) 


§/W-n’ 


and 


Pfi(xj)  Pfijxj)  _ Wl . _ p f2(xi) 


3=1 


fiXj) 


U I f(xi) 


= ^2wi3  - pn,  = 0 


so  that 


1 n 

p = - ]£ WH  ■ 


nT=  i 


From  the  second  equation  in  (3.2.4), 


2 = S?=l  XjWii 
‘ E"=l  Wfc 


(3.2.3) 


(3.2.4) 


(3.2.5) 


(3.2.6) 


i = 1,2  . 


(3.2.7) 
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Thus,  from  (3.2.6)  and  (3.2.7),  the  (z/  -f  l)th  iteration  forms  for  p and  of  (z  = 1, 2) 
are 


p<-+i)  = ix>i? 

Y'n  x2wW 

( 2)(,+1)  = ^=,1...;.,  V._  i = 12 


where  is  an  estimate  based  on  the  parameters  of  the  iteration  of  the  dis- 
tribution, i.e., 

M _ 


WH  = 


ft*  (®i) 


where  f[v\xj\i  = 1,2  and  f^v\x  j)  are  the  density  function  evaluated  at  Xj  using 
the  parameter  estimate,  p^v\  cr2^  and  . The  iterations  continue  until  some  pre- 
determined criteria  has  been  met  which  stops  the  procedure.  The  above  procedure 
has  been  discussed  by  Behboodian  (1970),  and  Hosmer  (1971,  1973a),  and  others. 
Let  the  above  maximum  likelihood  estimates  be  denoted  by  p and  <x»2,  i = 1,2. 


3.3  Modified  Maximum  Likelihood  Estimation 

A new  estimation  method  is  developed  for  more  accurate  estimates  of  a\  and 
proportion  p.  Since  the  main  objective  is  estimating  a2,  the  new  estimation  method 
uses  a truncated  distribution,  which  contains  more  information  about  the  smaller 
variance,  a2.  The  procedure  involves  three  steps  and  is  explained  below. 

In  the  first  step,  <7i  is  estimated,  say  a\,  by  using  doubly  truncated  normal 
distributions.  For  the  truncation  points,  a multiplication  of  the  standard  maximum 
likelihood  estimate  di  described  in  Section  3.2  is  used. 
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In  the  second  step,  estimates,  <r$  and  p*,  are  obtained  for  cr2  and  p respectively, 
based  on  a mixture  of  two  normal  distributions  with  given  &\  estimate  from  the 
first  step. 

In  the  third  step,  estimate,  <r*,  is  obtained  for  a i,  based  on  a truncated  mixture 
of  two  normal  distributions  with  given  and  p*  estimates  from  the  second  step. 
Truncation  points  axe  the  same  as  those  of  the  first  step.  Detailed  mathematical 
explanations  follow. 

3.3.1  Estimation  of  <j\  from  Doubly  Truncated  Single  Normal  Distribution 


The  density  function  for  a truncated  normal  random  variable  with  truncated 


points  at  ±xo  is 


/(*)  = 


(F(x0)  - F(-Xo))V^7T<7i 

i =4 


e aoi 


(2  F(x0)  — l)^i 

(2$(^)  _ 1)^7? 

where  F is  the  cumulative  distribution  of  N(0,<Ti),  $ is  the  cumulative  distribution 
of  N(0,  1),  and  (f>  is  the  probability  density  function  of  N(0,1).  The  log  likelihood 
function  is 


L = ln  II  f(xi)  = 12 ln  /(*j) 

j= i j= i 

= ln(2$(— ) — 1)  — ln  <xi  — ln  4>{—)  • 


j=x 


Differentiation  with  respect  to  gives 
dL 


dL  _ y.  f 1 Xj  ) 

dal  ~ ,ti\24(a)-l  2tr?  + la'  ] 
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2*?  12#(“)-1  W + v?  I 


since 


**(£)  _ 


d ,XQ\]  A(* On  _ x0  i/SQi 

5(7i  [5<Ti  CTX  2 (Tj  O'!  ’ 


Then  = 0 leads  to 

C7<rf 


i (ZU*V 


aS  = 


1-2^ 


n 


where,  | = ?,  0 = fffl  ■ 

<^i  2$(£)  - 1 


Thus,  the  ( v + l)th  iteration  form  for  estimation  of  o\  is 

(crHl/+1'> 1 ( ^i=i  xi  \ 

Kx)  ~ 1 -2q^)f^  n h 


where 


cH  - II. 

c ~ M ’ 
<t\ 


9 = 


<f>U{v)) 


(3.3.1) 


2$(£W)-1’ 

and  and  are  the  probability  density  function  and  cumulative  distri- 

bution of  N(0,  1)  at 


3.3.2  Estimation  of  a2  and  Proportion  p from  Mixture  of  Two  Normal  Distributions 
Given  ct\ 

The  density  function  of  a mixture  of  two  normal  distributions  with  a pro- 
portion of  p from  N(0,  cTj)  and  a proportion  of  (1  — p ) from  N(0,  o\ ) is 

/(»)  = * + 0^.3 

V/27T0'i  \f2Tra2 

= pfi{x)  + (l  - p)fl{x) 

where  0 < p < 1 and 


/<(*)  = 


\/2Tr<Ji 


e2ai  , * = 1,2. 
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The  log  likelihood  function  is 


L = ln  n /(*>)  = 2 ln  /(®i)  = £ ln  (p/i(*j)  + (i  - P)Mxj)) 

j=i  j= l j=i 

Differentiation  with  respect  to  the  proportion  p and  cTj  gives 


Also 


/i(gj)  ~ /a(aj) 

0P  ^ /(**) 


= A(l-p) 

^ /(®i)  \ 2v^rcr; 

= £(l-p)J  J 


ft  /(X,)  I 2^+4/jM 

1 v-'/  2 2\(1  — p)fl(xj) 

/(*,)  • 


ft  = 0 

dp 

ft  = o 

provide  right  hand  side  of  (3.3.2)  and  (3.3.3)  to  zero. 
Let  pi  = p and  p2  = 1 — p and  define 


Wij  = Pity**3}  * = 1.2,  j = 1, . . . ,n. 

f(xj) 

Then  the  second  equation  of  (3.3.4)  results  in 

r *,2  Ei=i*iwSi 
(<T2)  - • 

By  the  same  approach  as  in  equation  (3.2.5)  and  (3.2.6),  p*  = £"=1  *°ij7 
Thus,  the  (u  + l)th  iteration  forms  for  a\  and  p are 


.2\(H-1) 

2) 

T.U  “S’ 

p("+1)  = 

1 (»>) 

’ 

ni=i 

(3.3.2) 


(3.3.3) 


(3.3.4) 
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where  w ^ is  an  estimate  based  on  the  parameters  of  the  iteration  of  the  dis- 
tribution, i.e., 

(»)  _ Pi^fi^jXj) 
ij  ~ 

where  = 1,2  and  f^v\x  j ) are  the  density  function  evaluated  at  Xj  using 

the  parameter  estimate,  p^v\  cr^u\  and  dj. 


3.3.3  Estimation  of  ai  from  Doubly  Truncated  Mixture  of  Two  Normal  Distributions 
Given  (72  and  p 

The  density  function  is 


Q(x)  — F(x0)-F(-X  o)^(Z) 

= f{x0)-f(-Xo)(pMx)  + 0-  - p)Mx)) 


where 


A(*)  = 


/»(*)  = 

f(x)  = p/i(i)  + (1  - p)fi(x) 

F(x0)  = /!So(p/i(®)  + (1  - p)f3(x))dx 

= pife)  + (i-p)*(s) 

F(-x0)  = 1 - F(x0) 

=►  F(i.)  - F(-xc)  = 2F(x0)  - 1 
The  log  likelihood  function  is 


I = In  n Q(*i)  =-E  l-(2F(x„)  - 1)  + E In  /(*,') 

0= i i=i  i=i 
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Since 


d\nQ(x) 

dal 


FM]  = 


dal 


2(&FM)  (£/(*)) 

2 F(x0)  - 1 + f(x) 


PXc 


Krr)> 


2 al  v a i 


[/(*)]  = ^a(p/i(*)  + (l  -p)/a(*)) 
^ / P_ 


9<7i  \V57r<Ji 


e^r 


differentiation  with  respect  to  gives 


dL 

dal 


E 

i=l 


d\nQ(xj) 


dal 


= E 


£U  (2F(Xo)  - 1)  V 2<T1  CT1 
1 

2^? 


'’•#*)]  + 1 


f(xj)2al^X^  (a 


2np  (x°\  x(x°\  i V''  Ph(xj)  f_2  _2\ 

SxZFT  W + k ’ “ l) 


Let 


t = 


x0 


0-1 


Then  |4 

oat 


= 0 gives 


2nP  fUf\  i v*  P/i(xi)/T 2 2\  n 

553^^  + 3 S' TKT*’'  “ 1 ° " 0 


Let  pi  = p and  p2  = 1 — p and  define 


(3.3.5) 


Wij  = Pi 


Mxi) 

f(xi) 


z — 1,2  j — 1 , . . . , 


n. 
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Then  equation  (3.3.5)  gives 


*2  _ 
°\  = 


En  , w*  x2 


2sj=l  Wlj  2F*(xa)-l 

Thus,  the  (i/  + l)th  iteration  form  for  cr\  is 

Vn  nSv)r2 

( ✓t2'\(*/+1)  — 2^j=i  w\j  xj 

' 1 ’ y>n  _ 2npf(-)»(f(-)) 

wlj  2F<v)(x0)-\ 


(3.3.6) 


where  w ^ is  an  estimate  based  on  the  parameters  of  the  iteration  of  the  dis- 
tribution, i.e., 

* /<->(*>) 

where  = 1,2  and  f^v\xj)  are  the  density  function  evaluated  at  Xj  using 

the  parameter  estimate,  cr2^\  p*  and  (crj)2. 


3.4  Simulation  Results 

Without  loss  of  generality,  assume  a\  = 1 and  o\  = k 2,  i.e.,  mixture  of  two 
normal  distributions  with  a proportion  p observations  from  N(0,  1)  distribution  and 
the  remaining  proportion  1—  p observations  from  N(0,  k 2)  distribution  is  considered. 
The  criteria  used  to  terminate  the  iterations  are  if  10  < v < 1000  and  |.L(0(,'+1))  — 
L(0^)\  <1  x 10-4  then  9 = 0(|/+1),  otherwise  9 = $(1000),  where  6 is  the  parameter 
to  estimate  and  L(9)  is  the  log  likelihood  function.  It  was  found  1 x 10-4  is  the 
best  compromise  in  computation  time  and  final  value  of  L (Hosmer,  1971). 

When  the  sample  size  is  large,  both  maximum  likelihood  estimation  method  and 
the  new  estimation  provide  accurate  estimates.  Thus,  the  comparisons  are  made  for 
a relatively  small  sample  size,  50,  and  for  k = 3,4  and  5.  The  square  root  of  mean 
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squared  errors  (SQMSE)  of  estimates  of  the  parameters  oq,  p and  <72  under  5000 
simulations  for  different  k,  p and  truncation  points  are  given  in  Table  3.1  through 
Table  3.3.  The  table  values  show  the  new  estimation  method  gives  a smaller  mean 
squared  errors  for  a*  and  p*,  but  not  for  a*.  For  example,  the  SQMSE  of  new 
estimate  for  sample  size  50,  k = 3,  p = 0.975  and  4di  truncation  point  is  0.1753, 
while  the  corresponding  SQMSE  of  maximum  likelihood  estimate  <7X  is  0.2170  (Table 
3.1).  Under  the  condition  of  truncation  point  5<7X,  k = 5,  and  p — 0.99,  the  SQMSE 
of  new  estimate  a*  is  0.1678,  while  that  of  maximum  likelihood  estimate  b\  is  0.2135 
(Table  3.1).  The  SQMSE  of  proportion  p*  and  p for  k = 4,  p = 0.975  and  truncation 
point  4<7i  are  0.1402  and  0.2672  (Table  3.2).  The  SQMSE  of  new  estimate  crj  does 
not  provide  smaller  values  than  that  of  maximum  likelihood  estimate  <72  (Table 
3.3).  For  example,  for  k = 4,  p = 0.975  and  truncation  point  4d1,  the  SQMSE  of 

is  2.5454,  while  that  of  <72  gives  2.5308  (Table  3.3).  Since  the  main  objective  is 
accurate  estimation  of  parameters  <7i  and  p,  the  new  estimation  method  accomplish 
our  goal.  Simulation  results  also  show  that  there  are  no  big  differences  for  the 
truncation  point  with  4,  5,  and  6 times  <7i.  Considering  more  information  on  the 
main  body,  4 times  <7j  will  be  a good  choice  for  clinical  practice. 

3.5  An  Example 

To  demonstrate  the  above  estimation  method  in  a mixture  distribution,  we  use 
the  clinical  periodontal  data  to  examine  the  measurement  error  of  using  the  standard 
electronic  probe,  to  measure  the  clinical  probing  attachment  level  under  various 
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Table  3.1: 


Square  root  of  mean  squared  error  of  a*  when  sample  size=50  based  on 
5000  simulations 


Truncation  point  at 


k p 

3(7! 

4<Ji 

5(7! 

6(7i 

m.l.e. 

3 0.9 

0.2060 

0.2015 

0.1958 

0.1950 

0.2169 

0.925 

0.2062 

0.1863 

0.1780 

0.1716 

0.2173 

0.95 

0.2033 

0.1882 

0.1765 

0.1596 

0.2168 

0.975 

0.2000 

0.1753 

0.1648 

0.1737 

0.2170 

0.99 

0.1989 

0.1765 

0.1734 

0.1597 

0.2190 

4 0.9 

0.1787 

0.1676 

0.1681 

0.1697 

0.1869 

0.925 

0.1834 

0.1679 

0.1675 

0.1560 

0.1922 

0.95 

0.1840 

0.1788 

0.1648 

0.1518 

0.1988 

0.975 

0.1939 

0.1845 

0.1734 

0.1466 

0.2099 

0.99 

0.1975 

0.1936 

0.1533 

0.1622 

0.2174 

5 0.9 

0.1614 

0.1564 

0.1597 

0.1535 

0.1658 

0.925 

0.1635 

0.1521 

0.1553 

0.1504 

0.1718 

0.95 

0.1764 

0.1562 

0.1555 

0.1417 

0.1832 

0.975 

0.1845 

0.1684 

0.1664 

0.1538 

0.2003 

0.99 

0.1930 

0.1702 

0.1678 

0.1517 

0.2135 
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Table  3.2: 


Square  root  of  mean  squared  error  of  p*  when  sample  size=50  based  on 
5000  simulations 


Truncation  point  at 


k p 

3(7! 

4<7i 

5(7i 

6di 

m.l.e. 

3 0.9 

0.1394 

0.1040 

0.0836 

0.0727 

0.2197 

0.925 

0.1493 

0.1128 

0.0896 

0.0778 

0.2364 

0.95 

0.1630 

0.1252 

0.1030 

0.0925 

0.2573 

0.975 

0.1816 

0.1454 

0.1248 

0.1142 

0.2841 

0.99 

0.1968 

0.1633 

0.1446 

0.1352 

0.3051 

4 0.9 

0.1103 

0.0818 

0.0656 

0.0576 

0.1703 

0.925 

0.1245 

0.0939 

0.0743 

0.0636 

0.1933 

0.95 

0.1428 

0.1105 

0.0893 

0.0797 

0.2227 

0.975 

0.1720 

0.1402 

0.1201 

0.1101 

0.2672 

0.99 

0.1933 

0.1622 

0.1421 

0.1332 

0.2980 

5 0.9 

0.0893 

0.0692 

0.0577 

0.0520 

0.1327 

0.925 

0.1030 

0.0791 

0.0634 

0.0545 

0.1574 

0.95 

0.1248 

0.0986 

0.0807 

0.0709 

0.1940 

0.975 

0.1614 

0.1316 

0.1136 

0.1036 

0.2473 

0.99 

0.1877 

0.1576 

0.1386 

0.1302 

0.2899 
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Table  3.3: 


Square  root  of  mean  squared  error  of  a % when  sample  size=50  based  on 
5000  simulations 


Truncation  point  at 


k p 

3&i 

4<ti 

5&! 

6<Ti 

m.l.e 

3 0.9 

1.3827 

1.4198 

1.4308 

1.4536 

1.3674 

0.925 

1.5046 

1.5238 

1.5236 

1.5415 

1.4822 

0.95 

1.6260 

1.6330 

1.6302 

1.6370 

1.5991 

0.975 

1.7836 

1.7675 

1.7640 

1.7639 

1.7368 

0.99 

1.8967 

1.8709 

1.8688 

1.8637 

1.8385 

4 0.9 

1.7766 

1.8055 

1.8139 

1.8305 

1.7949 

0.925 

1.9821 

1.9966 

1.9931 

1.9976 

1.9957 

0.95 

2.2443 

2.2476 

2.2348 

2.2333 

2.2441 

0.975 

2.5533 

2.5454 

2.5321 

2.5283 

2.5308 

0.99 

2.7821 

2.7612 

2.7523 

2.7489 

2.7377 

5 0.9 

2.1284 

2.1443 

2.1507 

2.1623 

2.1674 

0.925 

2.4420 

2.4529 

2.4479 

2.4467 

2.4750 

0.95 

2.8343 

2.8329 

2.8212 

2.8116 

2.8504 

0.975 

3.3001 

3.2842 

3.2727 

3.2627 

3.2892 

0.99 

3.6788 

3.6557 

3.6452 

3.6395 

3.6413 
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conditions.  The  patients  were  classified  into  three  groups;  healthy,  gingivitis,  and 
periodontitis.  Ten  subjects  were  selected  from  each  group  and  36  typical  sites  per 
subject  were  carefully  probed  three  times  and  the  measurements  were  electronically 
recorded.  The  number  of  repeated  probings  was  limited  to  3 because  this  was 
considered  to  be  the  maximum  number  a patient  can  tolerate  during  one  visit. 
Details  of  this  experiment  are  described  in  Yang  et  al.  (1992). 

Since  the  true  attachment  level  at  each  site  is  unknown,  the  true  error  of  each 
probing  is  also  unknown.  But  the  error  distribution  can  be  approximated.  Let 
X\,X2  and  23  be  the  three  i.i.d.  random  variables  from  a mixture  of  two  normal 
distributions  with  a proportion  p from  iV(/z,  cr2)  and  q{=  1 — p)  proportion  from 
N(p,  k2a2).  Then  x>  — 2 j,  for  any  i,  j=l,  2,  3,  has  the  following  distribution. 

P2N( 0, 2cr2)  + 2pqN(0,  ( k 2 + l)a2)  + q2N( 0, 2/bV).  (3.5.1) 

Since  q is  small  (usually  less  than  0.05),  the  last  term  associated  with  q2  is  hardly 
observable.  Thus,  the  histogram  should  look  like  a mixture  of  two  normals  with  the  0 
mean  and  different  variances.  Figure  1.1  was  constructed  by  these  xl—Xj  differences. 
The  three  unknown  parameters  p,  <j2,  and  k are  estimated  by  the  previous  section’s 
modified  estimation  method  using  the  truncation  point  as  four  times  the  standard 
maximum  likelihood  estimate  di.  The  estimated  values  are  in  Table  3.4.  Table  3.5 
shows  the  simplification  of  Table  3.4.  The  second  through  fourth  columns  in  Figure 
3.1  illustrate  the  simulated  error  of  360  sites  by  the  estimated  error  distribution 
from  the  three  groups.  It  can  be  seen  that  p by  using  option-3  is  better  than  the 
mean  estimate  from  two  or  three  repeated  measurements.  The  first  column  figures 
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Table  3.4: 


Table  3.5: 


Modified  maximum  likelihood  estimates  of  parameters  in  (3.5.1)  for  3 
groups  of  dental  patients 


Group 

(y/2a)* 

(P2)* 

(^PTT a)* 

Healthy 

0.186 

0.811 

0.529 

Gingivitis 

0.182 

0.721 

0.563 

Periodontitis 

0.265 

0.875 

1.013 

Modified  maximum  likelihood  estimates  of  parameters  of  error  distri- 
bution for  3 groups  of  dental  patients 


Group 

<7* 

P 

k* 

Healthy 

0.13152 

0.9 

3.89 

Gingivitis 

0.12869 

0.849 

4.259 

Periodontitis 

0.18738 

0.9354 

5.313 
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(a)  x3-x(3)bar 


(a)  x2bar 


(a)  x3bar 


(a)  opt- 3 


(c)  x3-x(3)bar 


Figure  3.1:  Frequency  of  errors 

Note:  (a),  (b),  and  (c)  represent  healthy,  gingivitis,  and  periodontitis 
groups,  respectively.  x3-x(3)bar  represents  x3— (xi+a^)/^  based  on  real 
data.  x2bar,  x3bar  and  opt-3  represent  mean  of  two,  three  observation 
and  estimate  of  (1)  based  on  simulation,  respectively. 
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in  Figure  3.1  show  the  differences  x3  — (xi  + x3)/2,  the  best  we  can  get  for  the 
original  measurement  error.  The  patterns  are  very  similar  to  those  in  the  second 
and  third  column. 

3.6  Conclusion 

In  Chapter  2,  it  has  been  found  that  the  option-3  measurement  scheme  is  a 
function  of  two  variances  and  proportions  in  a mixture  of  two  error  distributions, 
especially  the  smaller  variance  which  is  crucial  to  the  threshold,  cr6,  determination. 
Thus,  the  main  concern  was  the  accuracy  of  the  smaller  variance  and  its  proportion. 

This  chapter  describes  a new  estimation  method  for  two  variances  and  propor- 
tions in  a mixture  of  two  normal  distributions.  It  is  found  that  the  new  estimation 
method  provides  a more  accurate  estimate  of  the  smaller  variance  and  its  proportion 
for  small  samples. 


CHAPTER  4 

DETECTION  OF  CHANGES  IN  LONGITUDINAL  DATA 
4.1  Introduction 

As  mentioned  in  Chapter  1,  detection  of  changes  in  longitudinal  data  is  often 
used  to  evaluate  the  progress  of  a disease.  For  example,  longitudinal  data  of  attach- 
ment level  (AL)  can  be  used  for  periodontal  disease  detection.  Three  methods  of  de- 
tecting AL  change  which  are  suggested  by  Haffajee  et  al.  (1983)  are  running  median, 
regression,  and  tolerance  methods.  Aeppli  and  Pihlstrom  (1989)  introduced  another 
method,  cusum,  and  compared  it  with  running  median  and  regression  method  under 
the  cases  of  the  gradual  or  burst  deterioration  model  in  disease.  They  claim  that 
the  cusum  and  running  median  methods  are  better  at  detecting  AL  change  than 
regression.  One  possible  explanation  for  poor  detection  of  change  by  the  regres- 
sion method  is  that  Aeppli  and  Pihlstrom  (1989)  do  not  use  the  same  assumption 
for  the  three  methods:  distribution  of  measurement  error  is  assumed  known  in  the 
running  median  and  cusum  methods,  but  not  for  regression.  Yang  et  al.  (1991) 
have  performed  a similar  comparison  under  the  same  assumption  that  distribution 
of  measurement  error  is  known  for  all  three  methods,  and  they  claim  that  regression 
is  better  in  most  cases.  There  is  another  discrepancy  in  the  assumed  error  structure 
in  the  two  papers.  Measurement  error  structure  is  a normal  distribution  in  Yang  et 
al.  (1991)  and  a logistic  distribution  in  Aeppli  and  Pihlstrom. 
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In  this  chapter,  the  regression,  running  median,  and  cusum  methods  are  recon- 
sidered and  compared  under  the  same  logistic  error  distribution  using  different  types 
of  models  for  change  in  deterioration.  It  is  found  that  the  regression  method  de- 
tects changes  more  quickly  than  the  running  median  and  cusum  methods  under  the 
logistic  error  distribution.  Four  new  methods  for  detecting  changes  are  introduced 
and  compared  with  the  regression  method  under  the  normal  error  structure.  After 
the  efficient  methods  for  detecting  changes  are  decided,  the  thresholds  for  given  a 
total  visit  number  and  specificity  are  provided.  Here,  specificity  is  defined  as  the 
probability  of  not  making  a false  alarm. 

4.2  Comparison  of  Regression,  Running  Median,  and  Cusum 

In  this  section,  the  regression,  running  median  and  cusum  methods  are  com- 
pared with  simulation  size  10,000.  Data  are  generated  based  on  logistic  random 
error  structure,  which  is  the  same  as  that  of  Aeppli  and  Pihlstrom  (1989).  Brief 
descriptions  of  data  structure  and  the  three  methods  are  given  here.  For  more  de- 
tailed explanations,  see  Haffajee  et  al.  (1983),  Aeppli  and  Pihlstrom  (1989),  and 
van  Dobben  de  Bruyn  (1963). 

4.2.1  Data  Structure 

Let  Yi,  Y2, . . . , Y„  be  the  longitudinal  data  over  time  t,  t = 1, 2, . . . , n.  Data  are 
created  under  the  same  model  as  that  of  Aeppli  and  Pihlstrom  (1989).  The  model 
is  Yt  = Xt  + Et  = Xa  -f  Dt  + Et,  where  X0  is  the  true  baseline  value,  Dt  is  the 
change  amount  at  visit  t,  and  Et  is  the  measurement  error.  True  baseline  value  is 
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estimated  by  the  average  of  two  values  which  are  randomly  selected  from  a uniform 
distribution  with  range  from  3 to  5 units.  In  this  section,  unit  measurement  is  mm. 
For  the  notation  convenience,  we  use  values  without  mm.  For  example,  5 means  5 
mm.  Change  amount  is  determined  by  the  types  of  change  models.  Three  possible 
deterioration  models  are  considered.  The  three  models  are 


Gradual  model  : 


Burst  model  : 


Yt  — X0  -f  Et,  t < t0 

Yt  = X0  + (3(t  — tQ)  + Et,  t>t0 

Yt  = X0  + Et,  t < t0 

i 

Yt  = XQ  + A + Et,  t > ta 


Random  Deterioration  model  : 


Yt  — X0  + Et,  t < tQ 

(4.2.1) 

Yt  = Xt~i  + 7t  + Et,  t >t0 


where  t0  is  the  change  point  and  (3  and  A are  determined  by  change  amount,  and  7* 
is  randomly  selected  from  uniform  distribution  such  that  It  equals  to  a given 


total  change  amount.  The  first  two  models  have  been  considered  in  most  dental 
studies  (Haffajee  et  al.,  1983;  Aeppli  and  Pihlstrom,  1989;  and  Yang  et  al.,  1991). 
The  third  model  is  added  because  we  feel  it  is  as  reasonable  as  the  previous  two. 
Measurement  errors  are  generated  from  a logistic  distribution  which  is  the  same 
error  generation  as  Aeppli  and  Pihlstrom’s.  The  variances  are  set  up  0.5  for  most 
3 unit  measurement  and  increase  linearly  by  0.075  per  unit  beyond  3;  for  example, 
if  the  measurement  is  4,  then  variance  is  assigned  0.575(0.5  + 0.075  = 0.575). 

The  three  methods  compared  in  Aeppli  and  Pihlstrom  (1989)  are  described  as 


follows. 
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4.2.2  Regression  Method 

The  regression  method  assumes  that  Yt  = a + /3t,  where  a is  the  intercept  and  (3 
is  the  slope.  The  slope(/3)  is  determined  by  the  least  squares  method  after  each  new 
observation  and  tested  against  zero  at  a preset  significance  level,  and  the  rejection 
of  the  null  hypothesis  of  zero  slope  is  considered  as  an  indication  of  change.  In 
mathematical  notation,  change  is  declared  if 


sd(/3) 


> Cl 


(4.2.2) 


where  sd  represents  standard  deviation,  Ci  is  the  threshold  value  determined  by 
preset  significance  level,  and  /?  is  the  least  squares  estimate  of  (3.  Aeppli  and 
Pihlstrom  (1989)  use  sd(/3)  as 

1/2 


° -<&]*/(" -2)} 


where  an  = 


n(n+  l)(n  - 1)  j1/2 


and  a is  the  least  squares  estimate  of  a.  In  the  latter  part  of  this  chapter,  sd  of 
measurement  error  is  assumed  known 


sd(/3)  = cr /a* 


(4.2.3) 


where  a is  the  standard  deviation  of  random  error. 


4.2.3  Running  Median  Method 

The  running  median  sequence  is  defined  as  the  sequence  of  medians  of  three 
successive  observations.  A change  is  claimed  if  the  absolute  value  of  the  difference 
between  the  median  of  the  last  three  observations  and  the  baseline  is  greater  than 
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a preset  amount.  In  mathematical  notation,  the  decision  can  be  represented  by 

declare  change  if  |median(Y^, h^-2)  — baseline|  > C2  (4.2.4) 


where  C2  is  a preset  value. 

4.2.4  Cusum  Method 

Differences  from  baseline  are  added  to  obtain  a sequence  of  sums.  If  the  sums  are 
greater  than  or  less  than  predetermined  values,  a change  is  claimed.  More  detailed 
explanation  is  in  Aeppli  and  Pihlstrom  (1989),  Brooks  and  Evans  (1972),  and  van 
Dobben  de  Bruyn  (1963).  A brief  mathematical  explanation  is  as  follows. 

Let 


So  — 0,  S0  — 0, 

S„  = max(0,  + Yn-  fi-k), 

S~  = min(0,  5“_j  + (Yn-  fi)  + k),  for  n > 0,  (4.2.5) 


then  change  is  declared  if  5+  > h or  S~  < —h, 

where  fi  is  the  baseline,  and  k and  h correspond  to  the  slope  of  V-mask  and  height 
of  V-mask,  respectively  (pp.  14,  Van  Dobben  de  Bruyn,  1963). 


4.2.5  Simulation  Result 

Random  numbers  are  generated  based  on  the  paper  by  Park  and  Miller  (1988), 
and  a comparison  of  simulation  results  is  based  on  size  10,000.  In  this  section,  mea- 
surement standard  deviation  a is  assumed  known;  thus,  without  loss  of  generality 
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a = 1 is  used  in  the  simulation  study.  It  is  assumed  that  the  data  are  collected  by 
one  measurement  at  each  visit.  So  the  ith  visit  measurement  corresponds  to  the  ith 
observation  in  the  data  set.  Specificity  and  sensitivity  defined  below  are  used  for 
the  comparison  of  above  three  methods. 

• Specificity  is  defined  as  the  probability  of  claiming  no  changes  when  there  are 
no  changes,  i.e.,  the  probability  of  not  signaling  a false  alarm.  Hence  1 — a is 
equal  to  specificity  if  a is  the  significance  level. 

• Sensitivity  is  the  probability  of  claiming  change  when  there  is  a change. 

Specificity  of  approximately  95%  over  the  22  visits  is  used  to  determine  the  preset 
critical  value,  i.e.,  threshold  value.  Under  this  condition,  Aeppli  and  Pihlstrom 
(1989)  found  that  c2  in  (4.2.4)  is  three  times  the  standard  deviation  (a)  for  the 
running  median  method,  ci  = 2.8  for  regression  in  (4.2.2),  and  k = 2.4a  and  h = 
2.5 a for  cusum  in  (4.2.5).  However,  our  simulations  indicated  that  when  using  the 
above  critical  values,  the  specificity  is  not  close  to  95%  . They  turned  out  to  be  96%, 
90%,  and  96.6%  for  running  median,  regression,  and  cusum  method,  respectively. 
Better  threshold  for  95%  specificity  over  22  visits  are  2.795a  for  running  median, 
3.115  for  the  regression  test,  and  2.23a  and  2.5a  for  k and  h,  respectively,  in  cusum 
test  (Figure  4.1). 

Under  these  threshold  values,  i.e.,  under  the  same  95%  specificity  value  over  22 
visits,  the  sensitivities  of  each  test  with  respect  to  a burst  model  of  2.0u  or  3.0a 
change  at  the  second  visit  show  in  Figure  4.2  and  Figure  4.3  that  the  regression 
test  identifies  changes  more  quickly  than  cusum  or  running  median  methods,  which 


56 


indicates  that  the  results  of  Aeppli  and  Pihlstrom  (1989)  are  biased  against  the 
regression  method  because  they  assumed  that  the  variance  a2  is  unknown  in  the 
regression  method  but  known  in  the  other  two  methods.  For  example,  by  the 
fourth  visit,  for  the  burst  change  of  2.0a,  regression  identifies  the  changes  46%  while 
running  median  identifies  39%,  and  cusum  36%  . By  the  tenth  visit  for  the  burst 
change  of  2.0a,  the  sensitivities  of  the  regression,  running  median,  and  cusum  are 
78%,  64%,  and  69%,  respectively.  The  sensitivities  of  regression,  running  median, 
and  cusum  by  the  fourth  visit  for  the  burst  change  of  3.0<r  are  89%,  84%,  and  84%, 
respectively.  By  the  tenth  visit  for  the  burst  change  of  3.0a,  regression,  running 
median,  and  cusum  detect  99%,  96%,  and  97%,  respectively. 

The  sensitivities  for  the  gradual  deterioration  over  time  with  (3  = 0.3a  and  0.5a 
per  visit  also  show  that  regression  test  is  more  sensitive  in  detecting  changes  than 
other  tests  (Figure  4.4  and  Figure  4.5),  which  agrees  with  Yang  et  al.  (1991).  For 
example,  by  the  sixth  visit  for  the  gradual  change  of  0.3a  (total  change  = 1.2a), 
regression  identifies  9%  of  the  changes,  while  running  median  identifies  6%  and 
cusum  5%.  The  sensitivities  of  regression,  running  median,  and  cusum  by  the  tenth 
visit  for  the  gradual  change  of  0.3a  (total  change  = 2.4a)  are  75%,  47%,  and  49%, 
respectively.  By  the  sixth  visit  for  the  gradual  change  of  0.5a  (total  change  = 2.0 
a),  regression  identifies  31%, while  running  median  identifies  19%  and  cusum  18%. 
The  sensitivities  of  regression,  running  median  and  cusum  by  the  tenth  visit  for  the 
gradual  change  of  0.5a  (total  change  = 4.0a)  are  99%,  95%,  and  96%,  respectively. 
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The  sensitivities  for  the  random  deterioration  model  show  similar  results  as  the 
gradual  change  model  (Figure  4.6).  By  the  tenth  visit,  for  random  deterioration 
change  with  total  change  4.0a  at  the  22nd  visit,  regression  detects  35%,  while  run- 
ning median  identifies  25%  and  cusum  20%.  Overall,  regression  identifies  changes 
more  quickly  than  the  other  two  methods,  running  median  and  cusum. 

4.3  New  Methods  in  Change  Detection 

It  has  been  demonstrated  by  Yang  et  al.  (1991)  that  using  cusum  is  inefficient 
when  the  baseline  attachment  level  has  to  be  estimated.  Thus,  there  is  no  reason 
to  use  cusum  method  when  the  baseline  is  unknown.  However,  the  running  median 
method,  which  uses  only  the  last  three  measurements  and  the  baseline,  can  be 
improved  by  incorporating  more  data  for  change  detection.  This  idea  is  to  identify 
the  change  point  first  and  then  use  the  difference  of  means  on  the  two  sides  of 
this  change  point  to  assess  the  possibility  of  change.  This  method,  referred  to 
later  as  single  change  with  change  point  detection  method,  and  its  counterpart  in 
regression  method  will  be  defined  and  investigated  in  this  section  (Hinkley,  1970; 
Quandt,  1958;  and  Zacks,  1983). 

4.3.1  Data  Structure 

Data  are  generated  under  the  same  model  of  section  4.2  with  normal  random 
error  instead  of  logistic  random  error.  Let  xi,  ij, . . . , xn  be  a longitudinal  data  set 


over  time.  The  data  set  is  created  as  follows. 
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Figure  4.1:  Specificities  of  three  methods 
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Figure  4.2:  Sensitivities  of  three  methods  for  burst  model  with 
change  point  at  visit  no.  2 and  2.0<7  change  amount 
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Figure  4.3:  Sensitivities  of  three  methods  for  burst  model  with 
change  point  at  visit  no.  2 and  3.0<r  change  amount 
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Figure  4.4:  Sensitivities  of  three  methods  for  gradual  model  with 
change  point  at  visit  no.  2 and  0.3a  change  amount 
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Figure  4.5:  Sensitivities  of  three  methods  for  gradual  model  with 
change  point  at  visit  no.  2 and  0.5cr  change  amount 
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Figure  4.6:  Sensitivities  of  three  methods  for  random  deterioration  model  with 
change  point  at  visit  no.  2 and  4.0<r  total  change  amount 
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• Null  hypothesis  model:  Since  there  is  no  change  over  time,  ®j,  ®2> . . . , xn  are 
generated  as  independent  and  identically  distributed  (i.i.d)  standard  normal 
random  variables. 

• Gradual  model:  For  a given  change  point  i0  and  change  rate  (3 , 


for  i < iQ  — 1 


/ 3(i  - i0 ) + e;,  for  i > ia 


where  et)  i = 1,. . . ,n  are  i.i.d.  standard  normal  variables. 


• Burst  model:  For  a given  change  point  i0  and  change  amount  A, 


Xi  — < 


for  i < ia  — 1 


A + tj,  for  i > ia 
where  e*,  i = 1, . . . ,n  are  i.i.d.  standard  normal  variables. 


• Random  deterioration  model:  X{  = £ + ej  i = 1, 2, . . . , n where  is  i.i.d. 

standard  normal  variable  and  Si  represents  random  magnitude  of  deteriora- 
tion. In  our  simulation  study,  Si  = 0 for  i < i0  and  Si  is  determined  by  uniform 
random  variable,  in  which  total  change  amount  is  given,  i.e.,  X^t>*0  = given 

total  change  amount. 


In  simulation,  random  numbers  are  again  generated  by  the  method  of  Park  and 
Miller  (1988),  and  normal  variables  are  generated  by  Box  and  Muller  transformation 
(Box  and  Muller,  1958). 
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4.3.2  Single  Change  with  Change  Point  Detection  (SCCPD)  Method 


Single  change  with  change  point  detection  (SCCPD)  method  is  designed  to  fit 
the  model 


/ii  + ti  if  i < ia 


Xi  = < 


+ e»  if  i > ia 

where  ^ and  /r2  are  the  true  means  for  i < ia  and  i > i0 , and  e/s  are  i.i.d.  random 
variables.  The  change  point  is  estimated,  say  i0,  by  the  one  which  gives  the  smallest 
fitting  error 


sse  , = Efe-fc)!+Efe-W! 

j<*  0 >>*0 


where 

/ti  = average  of  Xj,  j < i„ 

/t2  = average  of  Xj , j >ia  . 

The  SCCPD  method  detects  changes  by  using  the  standardized  difference  of  two 
means  separated  by  iQ.  If  the  standardized  difference  of  two  means  is  greater  than 
a predetermined  value,  change  is  declared.  More  precisely,  the  test  statistic  for 
SCCPD  with  the  estimated  iQ  is 


A2  ~ Ai 
sd(A2  - A 1) 

where  sd(A2_ Ai)  is  the  standard  deviation  of  A2— Ai-  i-e->  square  root  of  uar(s)/[r-^--f 
n_tx  +1].  If  equation  (4.3.1)  is  greater  than  predetermined  value,  then  change  is 
claimed. 


(4.3.1) 
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4.3.3  Regression  with  Change  Point  Detection  (RCPD)  Method 


Regression  with  change  point  detection  (RCPD)  method  is  designed  to  fit  the 
model 


ft  + ^ 


if  i <iQ 


X{  = < 


a + &{}  — ia ) + if  i > i0 

where  /z  is  the  true  mean  for  i < i0,  and  a and  /3  are  constant  and  slope  of  linear 
regression  model  with  starting  point  iQ , and  e^,  i = 1, . . . , n are  i.i.d.  normal  random 
variables.  The  change  point  is  estimated,  say  ia,  by  the  one  which  gives  the  smallest 
fitting  error 


SSE2  = £(*J  - - i„)f 

j<i  o J>*o 

A 

where  /t  = average  of  x3,  j < iQ  and  a and  (3  are  the  least  squares  estimates  for 
linear  regression  model  with  starting  point  i0.  RCPD  method  detects  changes  by 
the  test  of  zero  slope(/3)  in  linear  regression  model.  Rejection  of  null  hypothesis  of 
zero  slope  is  considered  as  an  indication  of  change. 

It  is  obvious  that  the  regression  model  is  good  for  the  gradual  model  and  the 
SCCPD  is  good  for  burst  model,  but  in  practice  the  model  for  real  change  is  un- 
known. The  following  two  methods  try  first  to  identify  the  model  from  the  data, 
and  then  perform  the  change  detection  test  for  the  selected  model. 


4.3.4  Mixed  Method  of  Type  I 

The  mixed  method  of  type  I is  the  mixture  of  the  regression  method  and  the 
SCCPD  method.  The  first  step  is  to  choose  the  more  appropriate  one  from  two 
methods,  regression  and  SCCPD  methods,  for  the  given  data  set.  Cross-validation 
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method  is  used  for  the  selection.  If  the  choice  is  the  regression  method,  the  test 
of  zero  slope  is  performed  to  detect  the  change.  If  the  choice  is  SCCPD,  then 
the  standardized  difference  of  two  means  (4.3.1)  is  used  for  the  test  of  change. 
Mathematical  details  are  given  as  the  following. 

Let  data  be  Si,  x2,  ■ • • > *n-  Consider  the  two  models 


Regression  X{  = a + /?  i + e. 


and 


Burst 


Xi  = 


/ii  + if  i < i0 
//2  + Ci  if  i > i0 


for  i = 1 ,n. 


Predicted  value  of  the  ith  observation  in  the  regression  model  with  the  ith  observation 

A A 

deleted  is  denoted  as  x ^ -f  /?(;),  where  d(i)  and  (3^  are  the  least  square 

estimates  of  constant  and  slope  of  linear  regression  from  the  data  without  the  ith 
observation.  In  SCCPD  method,  when  is  deleted,  estimated  change  point  iQ  is 
the  one  which  gives  the  smallest 


S (*j-£ i)2+  £ (xi~M2 

and  the  predicted  value  is 

{£i(0  if  * < *o 
fo{i)  if  * > t0. 

where 

A i(i)  = average  of  xjt  j < i01  j / i 
h(i)  = average  of  a j > iot  j ± i . 
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For  the  selection  of  one  method,  predicted  sum  of  squares  (PRESS)  of  each  model 
is  considered.  PRESS  is  defined  as  XXX»  — £(*))2-  The  meth°d  with  the  smaller 
PRESS  is  selected  for  the  change  detection  test. 

4.3.5  Mixed  Method  of  Type  II 

The  mixed  method  of  type  II  is  the  combination  of  the  SCCPD  method  and  the 
RCPD  method.  The  process  for  detecting  changes  is  the  same  as  that  of  type  I.  The 
first  process  is  to  choose  the  appropriate  method  based  on  PRESS,  and  once  the 
better  method  is  selected,  the  testing  procedure  for  detecting  changes  is  the  same 
as  the  selected  method. 

4.3.6  Simulation  Result 

In  the  simulation,  random  errors  are  generated  based  on  a normal  distribution 
with  mean  zero  and  variance  one.  The  data  are  created  based  on  25  measurements 
which  correspond  to  the  baseline  and  24  monthly  measurements  for  two  years.  Based 
on  100,000  simulations,  it  is  found  that  regression,  SCCPD  and  the  mixed  method 
of  type  I identify  changes  more  quickly  than  other  methods.  The  reason  the  RCPD 
and  mixed  method  of  type  II  failed  in  this  competition  is  because  they  are  quite 
sensitive  to  the  location  of  the  real  change  point.  However,  the  RCPD  method  could 
not  detect  the  change  point  with  great  accuracy. 

To  determine  the  critical  value  for  the  test,  a specificity  of  approximately  95% 
over  25  visits  is  used.  Under  this  consideration,  the  chosen  thresholds  are  2.62  for 
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the  regression  method  and  3.172  for  SCCPD  method.  For  the  mixed  method  of 
type  I,  the  thresholds  are  2.72  for  regression  and  3.27  for  SCCPD. 

Under  the  same  specificity,  the  sensitivities  in  change  detection  by  regression, 
SCCPD  and  mixed  method  of  type  I are  given  in  Figure  4.7  through  Figure  4.10 
for  different  change  models.  Results  for  the  burst  model  with  2.0<r  and  3.0a  change 
amounts  at  the  three  change  points,  5th,  10th,  and  15th  visit,  are  given  in  Figure 
4.7  and  Figure  4.8.  By  the  eighteenth  visit  for  the  burst  change  of  2.0a  with  change 
point  at  the  fifth  visit,  the  regression  method  identifies  85%  of  the  changes,  while 
SCCPD  method  identifies  79%  and  mixed  method  of  type  I 83%  (Figure  4.7).  The 
sensitivities  by  regression,  SCCPD,  and  mixed  methods  of  type  I by  the  eighteenth 
visit  for  the  burst  change  of  2.0cr  with  change  point  at  the  tenth  visit  are  91.5%, 
92.7%,  and  92.8%,  respectively  (Figure  4.7).  By  the  eighteenth  visit  for  the  burst 
model  with  3.0a  change  amount  and  fifteenth  visit  change  point,  the  sensitivities 
by  regression  method,  SCCPD  method,  and  mixed  method  of  type  I are  89%,  99%, 
and  98.6%,  respectively  (Figure  4.8). 

Sensitivities  in  change  detection  of  the  gradual  change  model  with  change  amounts 
0.2cr  and  0.3a  at  three  different  change  points,  5th,  10th,  and  15th  visit  are  given 
in  Figure  4.9  and  Figure  4.10.  For  the  earliest  change  point  among  the  three,  the 
regression  method  detects  changes  more  quickly  than  SCCPD  method  and  mixed 
method  of  type  I.  For  example,  by  the  fifteenth  visit  with  change  point  at  the 
fifth  visit  and  0.2a  change  rate  (total  change  amount  is  2.0a),  regression  detects 
changes  55%,  while  SCCPD  detects  50%  and  mixed  method  of  type  I 54%.  When 
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the  change  occurs  in  the  middle  of  visits,  the  tenth  visit,  with  0.2a  change  rate, 
regression,  SCCPD,  and  mixed  method  of  type  I identify  changes  59%,  63%,  and 
63%  at  the  twentieth  visit  (total  change  amount  is  2.0a).  For  the  change  point  at 
the  fifteenth  visit,  with  0.3cr  change  rate,  sensitivities  at  the  twenty-fifth  visit  are 
93.5%,  97.8%,  and  97.3%  for  regression,  SCCPD,  and  mixed  method,  respectively. 

Results  for  the  regression  method  are  better  than  for  the  other  two  methods 
for  the  early  change  point  but,  for  the  later  change  point,  the  SCCPD  method  is 
superior  to  the  other  two  methods.  The  mixed  method  of  type  I consistently  follows 
just  behind  the  best  one  among  the  three  methods. 

Figure  4.11  shows  sensitivities  of  the  regression  method  with  total  change  amount 
4.0cr  at  the  25th  visit  in  gradual  and  random  deterioration  models.  It  shows  that 
there  is  no  difference  between  the  random  deterioration  and  gradual  models  when 
the  regression  method  is  used.  The  possible  explanation  is  that  the  random  deteri- 
oration model  indicates  gradual  breakdown  over  visit,  though  not  at  the  same  rate 
at  every  visit.  Thus,  if  the  total  change  amount  is  the  same  in  the  random  and 
gradual  models,  almost  the  same  sensitivities  can  be  expected. 

4.3.7  Threshold  Determination 

Once  the  most  efficient  methods  for  change  detection  are  discovered,  the  thresh- 
olds for  claiming  changes  still  may  not  be  simple.  Thus,  in  this  section,  thresholds 
for  three  selected  best  methods  are  provided.  In  most  cases,  therapeutic  action  is 
required  as  soon  as  there  is  enough  evidence  that  deterioration  has  occurred  at  one 
site.  Hence  the  actual  clinical  diagnosis  is  a multiple  sequential  decision. 
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Figure  4.7:  Sensitivity  of  burst  model  with  2.0a  change  amount 
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Figure  4.8:  Sensitivities  of  burst  model  with  3.0cr  change  amount 
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Figure  4.9:  Sensitivities  of  gradual  model  with  0.2<j  change  rate 
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Figure  4.10:  Sensitivities  of  gradual  model  with  0.3cr  change  rate 
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Figure  4.11:  Sensitivities  of  regression  method  with  4.0 a total  change  amount 
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Consequently,  the  total  number  of  visits  and  number  of  sites  examined  are  essential 
parts  in  this  decision. 

1 

Let  N denote  the  total  number  of  visits.  For  the  general  application  of  the  above 
three  methods,  regression,  SCCPD,  and  mixed  method  of  type  I,  the  thresholds  for 
N=7,  13,  25  and  different  specificity  levels  (=  1 — a if  a is  significance  level)  are 
given  in  Figure  4.12  and  Figure  4.13.  Then  N=7,  13,  25  correspond  to  a one  time 
measurement  for  initial  visit  plus  monthly  measurements  of  a half-,  one-,  and  two- 
year  programs.  The  x-axis  of  Figure  4.12  and  Figure  4.13  corresponds  to  log10a. 
Thresholds  for  other  value  of  N can  be  obtained  by  interpolation.  The  following 
formula  provides  a good  approximation  for  the  figures. 

Regression  : c=  -0.2216  + 0.0158N  + 2.13466  y/\log10a\  R2  = 0.9976 

SCCPD  : c = 0.13618  + 0.0265N  + 2.03215  ^ \log10ct\  R2  = 0.9913 

Mixed(type  I)  : Regression  : c=— 0.1181  + 0.01726N+2. 11074  \J\logwa\  R2  — 9969 
SCCPD  : c=0.2667  + 0.02813N+1.9862  J\log10a\  R2  = 0.9922 

where  R2  is  the  coefficient  of  multiple  determination  in  regression  analysis.  Thus, 
for  any  given  a and  N,  the  thresholds  can  be  computed  by  the  above  formula. 
For  example,  the  regression  method  with  N=20  and  a = 0.001  gives  an  x-axis 
value  of  log10  0.001  = —3  and  threshold  3.7917.  For  SCCPD,  the  threshold  for 
N=17  and  a = 0.005  gives  3.669  with  x-axis  value  log10  0.005  = —2.301.  The 
threshold  for  the  mixed  method  of  type  I for  N=22  and  a = 0.001  are  3.9175  and 
4.32576  for  the  regression  method  and  SCCPD  method,  respectively.  Once  the 
threshold  is  obtained,  a test  of  change  detection  can  be  performed.  For  example, 
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if  the  regression  method  is  used  for  N=20  and  a = 0.001,  then  change  is  declared 
whenever  /3/sd(/3)  > 3.7917. 

The  above  formula  and  Figure  4.12  and  Figure  4.13  also  can  be  used  for  multiple 
site  measurements.  Suppose  there  are  s sites  and  the  possible  total  visits  is  N.  The 
usual  pattern  is  that  change  is  declared  if  change  is  detected  in  at  least  one  site. 
For  such  multiple  site  measurements,  thresholds  can  be  determined  for  the  above 
formula  by  substituting  a by  1 — (1  — a)1/'  if  the  overall  significance  level  is  a.  Thus, 
the  x-axis  is  changed  to  log10(l  — (1  — a)1/*).  For  example,  the  regression  method 
with  N=22,  s = 60  and  over  all  a = 0.05  yields  an  x-axis  value  of  log10(l— 0.951/60)  = 
—3.068  and  threshold  3.865. 

Thresholds  for  a one  time  decision  at  the  end  of  patient  visits  are  much  easier 
to  determine  than  those  of  a sequential  decision,  since  the  test  statistic  value  needs 
to  be  computed  only  at  the  end  of  the  data.  In  the  regression  method  case,  the 
threshold  is  $ where  $ 1 is  the  inverse  of  a standard  normal  distribution.  In 

the  multiple  site  case,  the  threshold  is  4>-1(l  — (1  — a)1/*)  when  s is  the  number  of 
total  sites. 

4.3.8  Unknown  Variance 

Specificity  and  sensitivity  are  affected  by  an  under-  or  over-estimated  variance. 
In  particular,  sensitivity  is  reduced  by  an  over-estimated  variance.  Thus,  when 
it  is  difficult  to  measure  some  sites  due  to  strange  tooth  structure,  it  becomes 
unreasonable  to  use  the  same  variance  for  all  sites.  One  possible  way  to  resolve 
this  problem  is  to  use  a method  that  requires  no  assumption  on  the  knowledge 
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Figure  4.12:  Thresholds  for  regression  and  SCCPD 
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Figure  4.13:  Thresholds  for  mixed  method  of  type  I 
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of  variance(<r2).  When  the  regression  method  with  unknown  a is  used  in  a model 
other  than  the  gradual  model,  e.g.  the  burst  change  model,  confounding  of  lack  of  fit 
error  and  measurement  error  cause  severe  problems  in  sensitivity  as  demonstrated 
in  Aeppli  and  Pihlstrom  (1989).  We  suggest  two  measurements  per  site  at  each 
visit  to  be  used  to  reduce  the  loss  of  sensitivity.  In  this  case,  a can  be  estimated 
independently  from  the  lack  of  fit  error.  Let  xH  and  x,2  be  the  two  measurements  at 
ith  time.  Then  the  ith  time  measurement  for  regression  is  determined  by  the  mean 

+ xi:i 

2 


and  a model  independent  estimated  standard  deviation  is 

In  this  case,  the  test  statistic  /3/sd(/3 ) has  a t distribution  with  n degrees  of  freedom. 
Thus,  the  thresholds  should  be  a function  of  size  n and  decrease  with  the  rate  of 
the  t-table  value.  Let  the  threshold  value  for  this  case  be  £(n)c  where  £(n)  is  a 
ratio  of  t quantile  and  normal  quantile.  For  the  multiple  site  case,  a small  a value 
needs  to  be  considered.  Thus,  the  t quantile  is  in  the  neighborhood  of  a = 0.0001 
level.  Table  4.1  shows  £(n)  value.  The  other  part  of  the  threshold,  c,  is  determined 
by  simulation  for  given  a and  number  of  total  visits.  Results  for  N=7,  13,  25  are 
plotted  on  Figure  4.14  and  other  thresholds  can  be  obtained  by  interpolation.  The 
following  formula  provides  a good  approximation  for  the  figures. 

c = -2.93558  - 0.12848N  + 1.09387^+2.323565  - 0.922065£2+ 


0.046435N5  - 0.012936N52  + 0.19807953  R2  = 0.999 
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where  B = \logioa\  and  R 2 is  the  coefficient  of  multiple  determination  in  regression 
analysis. 

The  comparison  of  sensitivities  of  the  known  and  unknown  variance  cases  under 
the  same  95%  specificity  is  given  in  Figure  4.15  and  Figure  4.16.  Sensitivities  for 
the  burst  change  model  with  2. Oct  change  amount  axe  on  Figure  4.15.  The  figure 
shows  that  there  is  substantial  difference  in  sensitivities  if  there  is  an  early  change 
point.  But  for  a later  change,  the  difference  in  sensitivities  is  small.  For  example, 
the  sensitivities  by  the  tenth  visit  with  a fifth  visit  change  point  are  77.95%  and 
94.38%  for  unknown  and  known  variance,  respectively,  while  the  sensitivities  by  the 
twentieth  visit  with  fifteenth  visit  change  point  are  96.92%  and  98.05%  for  unknown 
and  known  variance  case,  respectively. 

Sensitivities  of  the  gradual  change  model  (Figure  4.16)  show  results  similar  to 
the  burst  change  model  (Figure  4.15).  The  loss  of  sensitivity  by  unknown  variance 
under  the  gradual  change  model  is  less  than  that  of  the  burst  change  model.  By 
the  fifteenth  visit  with  the  change  point  at  the  fifth  visit  and  0.2ct  change  rate  (2. Oct 
total  change  amount),  the  sensitivities  are  81.93%  and  88.08%  for  unknown  and 
known  variance  case,  respectively.  The  known  variance  case  detects  90%  of  the 
changes  by  the  25th  visit,  while  the  unknown  variance  case  identifies  89.91%  of  the 
changes  when  the  change  occurs  at  the  fifteenth  visit  (2. 0<t  total  change  amount). 
Though  there  is  a loss  of  sensitivity  in  the  estimated  variance  case,  the  loss  is  small, 
especially  for  the  later  change  point. 
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Table  4.1:  £(n)  value 


n 

t 

2 

19.00851 

3 

5.969889 

4 

3.504391 

5 

2.602048 

6 

2.157670 

7 

1.899191 

8 

1.732108 

9 

1.616083 

10 

1.531025 

11 

1.466204 

12 

1.415250 

13 

1.374190 

14 

1.340426 

15 

1.312187 

16 

1.288231 

17 

1.267659 

18 

1.249807 

19 

1.234170 

20 

1.220364 

21 

1.208085 

22 

1.197096 

23 

1.187204 

24 

1.178253 

25 

1.170115 

Threshold 


83 


Figure  4.14:  Thresholds  for  regression  method  under  unknown  a 
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Figure  4.15:  Sensitivities  for  regression  method  for  known  and  unknown  a in  burst 
change  model  with  2.0a  change  amount 
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Figure  4.16:  Sensitivities  for  regression  method  for  known  and  unknown  a in  gradual 
change  model  with  0.2cr  change  rate 
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4.4  Conclusion 

Many  statistical  methods  for  detecting  changes  in  longitudinal  data  analysis  have 
been  studied.  In  periodontal  research,  three  methods,  regression,  running  median, 
and  cusum  have  been  suggested.  Among  the  three  methods,  the  regression  method 
identifies  changes  more  quickly  than  running  median  or  cusum  method  under  the 
three  possible  change  patterns:  gradual,  burst,  and  random  deterioration. 

Four  new  statistical  methods  for  detecting  changes  in  longitudinal  data  are  in- 
troduced and  compared  with  the  regression  method.  It  is  found  that  the  regression 
method,  SCCPD  method,  and  the  mixed  method  of  type  I,  are  more  sensitive  in 
detecting  changes  than  other  methods. 

For  the  purpose  of  application,  thresholds  for  regression  method,  SCCPD  method, 
and  mixed  method  of  type  I are  provided  for  sequential  multiple  decisions  with  vari- 
ous numbers  of  visits  and  sites.  Also,  for  the  unknown  error  variance  case,  thresholds 
are  provided  when  two  measures  are  allowed  per  visit. 


CHAPTER  5 

DETECTION  OF  LONGITUDINAL  CHANGE  UNDER  OPTION-3 

5.1  Introduction 

In  Chapter  2,  the  option-3  scheme  is  developed  to  reduce  the  influence  of  outliers 
so  that  correct  and  early  detection  of  changes  can  be  made  in  longitudinal  data.  In 
Chapter  4,  a numbers  of  change  detection  methods  for  longitudinal  attachment  level 
data  are  discussed.  Thus,  in  this  section,  the  option-3  scheme  is  used  for  change 
detection  of  longitudinal  data  and  compared  with  fixed  two  replicated  measurement 
scheme  by  specificity  and  sensitivity. 

5.2  Comparison  of  Option-3  Scheme  and  Fixed  Two-Sample  Scheme 

It  has  been  found  that  three  methods,  regression,  SCCPD,  and  mixed  method 
of  type  I,  axe  the  best  among  a number  of  other  methods  suggested  in  Chapter  4. 
Thus,  comparisons  here  are  made  only  for  these  three  selected  methods. 

Let  Xi,x2,X3  be,  respectively,  the  first,  second,  and  third  observations  in  each 
visit  which  come  from  an  i.i.d.  mixture  of  two  normal  populations  with  proportion  p 
from  N(0,  a2)  and  q(=l-p)  proportion  from  N(0,  k2cr 2).  Without  loss  of  generality, 
assume  a = 1.  The  Fixed  two-sample  scheme  represents  two  measurements  taken  at 
each  visit  and  regards  the  average  of  these  two  measurements  as  the  estimate  of  that 
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visit,  i.e.,  (xi  + x2)/2.0.  In  this  case,  the  variance  of  the  fixed  two-sample  scheme  is 
(p  + qk2)/ 2.0.  The  option-3  scheme  follows  the  same  procedures  as  Chapter  2 and 
the  variance  is  given  in  Table  2.1  for  the  optimal  6. 

In  this  section,  p = 0.95,  k = 4 and  p = 0.975,  k = 5 are  selected  to  compare 
the  option-3  scheme  and  the  fixed  two-sample  scheme.  All  simulations  are  based 
on  100,000  simulated  samples. 

5.2.1  p — 0.95  and  k = 4.0  case 

For  p = 0.95  and  k = 4.0,  the  thresholds  for  95%  specificity  are  2.907  for 
the  regression  method  in  fixed  two-sample  scheme,  2.674  for  regression  method  in 
option-3  scheme,  4.38  for  SCCPD  method  in  fixed  two-sample  scheme,  and  3.473 
for  SCCPD  method  in  option-3  scheme.  The  thresholds  for  mixed  method  of  type 
I are  2.902  and  4.5  for  regression  and  SCCPD  in  fixed  two-sample  scheme  and 
2.75  and  3.57  for  regression  and  SCCPD  in  option-3  scheme.  Figure  5.1  shows  the 
specificity  of  all  methods.  Under  the  same  specificity  at  25th  visit,  sensitivities  are 
given  in  Figure  5.2  through  Figure  5.5. 

Figure  5.2  shows  the  sensitivities  of  the  gradual  model  with  change  rate  0.2a 
for  three  methods.  In  the  regression  method,  by  the  15th  visit  with  the  change 
point(total  change  amount  = 2.0<r)  at  the  fifth  visit,  option-3  detects  77.8%  of 
changes  while  the  fixed  two-sample  identifies  50%.  Also,  the  sensitivities  by  the 
25th  visit  with  change  point  at  the  15th  visit  are  80%  for  the  option-3  scheme  and 
54%  for  the  fixed  two-sample  scheme.  In  the  SCCPD  method,  sensitivities  for  the 
option-3  scheme  are  much  better  than  for  the  fixed  two-sample  scheme.  By  the 
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15th  visit  with  change  point  at  the  fifth  visit,  the  option-3  scheme  detects  60%  of 
changes  and  the  fixed  two-sample  scheme  identifies  12%.  The  sensitivities  by  the 
20th  visit  with  change  point  at  the  10th  visit  are  76%  for  option-3  and  20%  for  the 
fixed  two-sample  scheme.  In  mixed  method  of  type  I,  the  option-3  identifies  73% 
of  changes  by  the  15th  visit  with  the  change  point  at  the  fifth  visit,  while  the  fixed 
two-sample  scheme  detects  40%.  The  sensitivities  by  the  20th  visit  with  the  10th 
change  point  are  81%  for  the  option-3  and  42%  for  the  fixed  two-sample  scheme. 
Overall,  Figure  5.2  shows  the  option-3  scheme  identifies  changes  more  quickly  than 
the  fixed  two-sample  scheme  for  all  three  methods,  regression,  SCCPD,  and  mixed 
method  of  type  I. 

Figure  5.3  shows  the  sensitivities  for  the  burst  model  with  2.0cr  change  amount. 
In  the  regression  method,  the  option-3  scheme  identifies  92%  of  changes,  while  the 
fixed  two-sample  scheme  detects  70.8%  of  changes  by  the  12th  visit  with  change 
point  at  the  fifth  visit.  By  the  16th  visit  with  a tenth  visit  change  point,  the 
sensitivities  are  96.8%  for  option-3  and  83.7%  for  fixed  two-sample  scheme.  For 
the  SCCPD  method,  the  option-3  scheme  detects  change  much  earlier  than  the 
fixed  two-sample  scheme,  the  same  phenomenon  as  observed  for  the  gradual  change 
model.  For  example,  the  sensitivities  are  94.6%  for  option-3  and  44%  for  fixed  two- 
sample  scheme  by  the  15th  visit  with  change  point  actually  occurring  at  the  10th 
visit.  By  the  20th  visit  with  the  15th  visit  change  point,  option-3  detects  97.7% 
of  changes  and  the  fixed  two-sample  scheme  60%  of  changes.  In  mixed  method  of 
type  I,  option-3  and  fixed  two-sample  scheme  detect  88.7%  and  52.6%  of  changes 
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by  the  12th  visit  with  change  point  at  the  fifth  visit.  The  sensitivities  are  99%  for 
option-3  and  83%  for  fixed  two-sample  scheme  by  the  20th  visit  with  change  point 
at  the  tenth  visit.  Thus,  Figure  5.3  shows  the  option-3  method  identifies  change 
more  quickly  than  the  fixed  two-sample  scheme  when  the  change  is  represented  by 
a burst  model. 

Figure  5.2  and  Figure  5.3  show  the  SCCPD  method  identifies  changes  later  than 
other  methods.  However,  the  SCCPD  method  is  comparable  with  other  methods 
in  Figure  4.7  and  Figure  4.9.  The  possible  reason  for  their  difference  is  that  in 
Chapter  4,  change  detection  was  performed  under  the  data  without  outliers,  i.e., 
under  a single  normal  distribution.  However,  here  the  detection  was  performed 
under  the  mixture  of  two  normal  distributions  which  include  outliers.  Thus,  these 
results  indicate  that  the  SCCPD  method  is  quite  affected  by  the  existence  of  outliers 
and  does  not  perform  as  well  in  their  present. 

Since  the  option-3  scheme  identifies  changes  more  quickly  than  the  fixed  two- 
sample  scheme  for  all  three  methods,  we  need  to  compare  the  three  methods,  re- 
gression, SCCPD,  and  mixed  method  of  type  I,  under  the  option-3  scheme.  Figure 
5.4  and  Figure  5.5  show  the  sensitivities  of  three  methods  for  the  gradual  and  burst 
models.  The  sensitivities  of  the  three  methods  under  the  option-3  scheme  show 
that  the  three  methods  are  comparable.  The  regression  method  identifies  changes 
more  quickly  for  the  model  with  an  early  stage  change  point,  while  the  SCCPD 
and  mixed  methods  detect  changes  earlier  for  the  model  with  a later  change  point 
(Figure  5.4  and  Figure  5.5).  For  example,  by  the  15th  visit  with  change  point  at  the 
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fifth  visit,  regression,  SCCPD,  and  mixed  method  of  type  I identify  77.8%,  60.4%, 
and  73%  of  changes,  respectively,  in  the  gradual  change  model  with  0.2  a change 
rate.  But  the  sensitivities  of  regression,  SCCPD,  and  mixed  method  of  type  I are 
45%,  49.3%,  and  50.9%,  respectively,  by  the  23rd  visit  with  change  point  at  the  15th 
visit  for  the  gradual  change  model  (Figure  5.4).  Similarly,  in  the  burst  model  with 
2.0cr  change,  by  the  15th  visit  with  the  change  point  at  the  fifth  visit,  regression, 
SCCPD,  and  mixed  method  of  type  I identify  95%,  88.7%,  and  92.8%  of  changes, 
respectively.  By  the  22rd  visit  with  change  point  at  the  15th  visit,  the  sensitivities 
are  98.7%  for  regression,  99.4%  for  SCCPD,  and  99.4%  for  mixed  method  of  type  I 
for  the  burst  change  model  (Figure  5.5). 

5.2.2  p = 0.975  and  k — 5.0  case 

For  p = 0.975  and  k = 5.0,  the  thresholds  for  95%  specificity  at  the  25th 
visit  are  2.954  for  the  regression  method  in  the  fixed  two-sample  scheme,  2.646  for 
the  regression  method  in  the  option-3  scheme,  4.754  for  the  SCCPD  method  in 
the  fixed  two-sample  scheme,  and  3.333  for  the  SCCPD  method  in  the  option-3 
scheme.  The  thresholds  for  the  mixed  method  of  type  I are  2.88  and  4.954  for 
regression  and  SCCPD  in  the  fixed  two-sample  scheme,  respectively,  and  2.75  and 
3.423  for  regression  and  SCCPD  in  the  option-3  scheme,  respectively.  Figures  5.6 
through  5.9  show  specificities  of  the  three  methods  under  gradual  and  burst  change 
models  with  the  same  specificity  95%  at  the  25th  visit. 

Figures  5.6  through  5.9  show  the  same  patterns  as  Figures  5.2  through  5.5.  For 
example,  the  option-3  scheme  detects  changes  earlier  than  the  fixed  two-sample 
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scheme  (Figure  5.6  and  Figure  5.7).  Also,  sensitivities  for  the  three  methods  under 
the  option-3  scheme  gives  similar  results  as  previous  ones,  i.e.,  for  the  change  model 
with  an  early  change  point,  the  regression  method  identifies  changes  a little  earlier 
than  the  other  two  methods,  and  for  the  later  change  point,  the  regression  method 
is  shown  to  identify  change  a little  later. 

5.3  Conclusion 

The  option-3  and  fixed  two-sample  schemes  are  used  for  change  detection  in  lon- 
gitudinal data.  A mixture  of  two  normal  distributions  is  used  to  explain  outliers. 
Two  choices  of  mixed  combination,  p = 0.95,  A:  = 4 and  p = 0.975,  k = 5,  are  used 
for  the  simulation.  The  results  of  simulation  for  the  two  choices  have  demonstrated 
that  the  option-3  scheme  detects  changes  earlier  than  the  fixed  two-sample  scheme 
and  the  three  methods,  regression,  SCCPD,  and  mixed  method  of  type  I,  are  com- 
parable change  detection  method  even  under  data  with  outliers  and  the  option-3 
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Figure  5.1: 


Specificities  of  the  three  methods  under  fixed  two  sample  and  option-3 
schemes  with  p = 0.95  and  k = 4.0 
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Figure  5.2:  Sensitivities  of  three  methods  for  the  gradual  model  with 

0.2cr  change  rate  under  the  fixed  two-sample  and  option-3  schemes 
for  p = 0.95  and  k = 4.0 
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Figure  5.3:  Sensitivities  of  three  methods  for  the  burst  model  with 

2.0 a change  amount  under  the  fixed  two-sample  and  option-3  schemes 
for  p = 0.95  and  k = 4.0 
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Figure  5.4:  Sensitivities  of  three  methods  under  the  option-3  scheme  for  the 
gradual  model*with  0.2 o change  rate  for  p = 0.95  and  k = 4.0 
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Figure  5.5:  Sensitivities  of  three  methods  under  the  option— 3 scheme  for  the 
burst  model  with  2.0cr  change  amount  with  p = 0.95  and  k = 4.0 
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Figure  5.6:  Sensitivities  of  three  methods  for  the  gradual  model  with 

0.2(7  change  rate  under  the  fixed  two-sample  and  option-3  schemes 
for  p = 0.975  and  k = 5 
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Figure  5.7:  Sensitivities  of  three  methods  for  the  burst  model  with 

2.0 a change  amount  under  the  fixed  two-sample  and  option-3  schemes 
for  p = 0.975  and  k = 5 
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Figure  5.8:  Sensitivities  of  three  methods  under  the  option-3  scheme  for 
gradual  model  with  0.2 a change  rate  for  p = 0.975  and  k = 5 
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Figure  5.9:  Sensitivities  of  three  methods  under  the  option-3  scheme  for  the 
burst  model  with  2.0cr  change  amount  for  p = 0.975  and  k = 5 


CHAPTER  6 

CONCLUSION  AND  FUTURE  RESEARCH 


Change  detection  in  longitudinal  data  is  an  important  diagnostic  tool  in  the 
medical  sciences.  Changes  in  clinical  probing  attachment  level  over  time  are  used 
for  determining  the  progress  of  periodontal  disease  in  dentistry.  However,  outliers 
in  measurement  errors  can  have  a detrimental  effect  in  change  detection.  A mixture 
of  two  normal  distributions,  in  which  one  has  a usual  variation  and  the  other  a 
relatively  large  variation,  is  considered  as  the  model  for  this  type  of  data.  To  account 
for  the  influence  of  outliers,  an  option-3  scheme  is  proposed.  The  form  of  option-3 
scheme  is  given  in  equation  (2.2.1),  which  involves  a threshold  08  determined  by 
minimum  mean  squared  error.  Due  to  the  complexity  of  the  computation  for  the 
mean  squared  error,  a numerical  solution  for  8 is  obtained  and  the  result  is  given  in 
Figure  2.1.  Under  optimal  8 which  gives  minimum  mean  squared  error,  the  option-3 
method  has  advantages  over  a fixed  repeated  measurements  scheme  in  sample  size 
and  kurtosis.  Since  outliers  are  relatively  scarce,  outlier  distribution  is  difficult  to 
determine.  Robustness  of  the  option-3  scheme  against  long-tail  and  short-tail  error 
distributions  are  examined,  and  the  result  provides  the  property  of  robustness  for 
outlier  distribution. 

To  determine  the  optimal  threshold,  the  parameters  of  the  error  distribution, 
especially  the  variance  of  the  normal  data  and  proportion,  are  needed.  A new 
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estimation  method  for  parameters  of  mixture  of  two  normal  distributions  is  provided 
in  section  3.3.  The  new  estimates  of  variance  of  the  normal  data  and  proportion 
under  simulation  are  shown  in  Table  3.1  and  Table  3.2.  These  are  more  accurate 
than  the  standard  maximum  likelihood  estimates  when  the  sample  size  is  small. 

Currently  used  change  detection  methods  in  periodontal  disease  are  the  regres- 
sion, running  median,  and  cusum.  These  methods  were  compared  and  it  was  found 
that  the  regression  method  identifies  changes  more  quickly  than  other  methods. 
Four  new  methods  for  detecting  changes  were  developed  and  compared  with  regres- 
sion. The  result  showed  that  regression,  single  change  with  change  point  detection, 
and  mixture  of  two  methods  were  more  sensitive  than  other  methods  in  detection 
of  changes.  The  option-3  and  a fixed  two  repeated  measurements  scheme  were 
compared  under  the  selected  three  efficient  change  detection  methods.  The  results 
showed  that  the  change  detection  method  with  the  option-3  measuring  scheme  could 
identify  changes  more  quickly  than  the  one  with  fixed  two  repeated  measurements 
scheme. 

The  option-3  scheme  was  developed  based  on  a mixture  of  two  normal  distri- 
butions and  showed  a robust  scheme  against  nonnormal  outlier  error  distribution 
with  normal  distribution  in  the  main  body.  Sometimes,  the  main  body  of  the  dis- 
tribution itself  is  not  normally  distributed.  Thus,  for  future  research,  the  option-3 
scheme  can  be  extended  to  the  non-normal  main  error  distribution.  Also,  we  may 
consider  a scheme  for  general  sample  size,  not  just  for  two  or  three  samples.  In 
parameter  estimation,  the  new  estimation  method  was  quite  complicated;  also,  it 
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used  standard  maximum  likelihood  estimates  for  the  truncation  point.  A parameter 
estimation  method  which  requires  less  complex  procedures,  uses  a simple  truncation 
point,  and  gives  more  accurate  estimates  for  both  variances  is  another  interesting 
research  topic  in  the  future. 
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